Description Usage Arguments Details Value Note Author(s) References See Also Examples
Given the loglikelihoods for a set of models arranged in ascending order of size, the best model is selected using the BICq criterion for a specified size.
1 |
logL |
vector of loglikelihoods |
n |
sample size |
level |
probability |
mSize |
model sizes |
mComplex |
a complexity function |
See reference
khat |
dataframe with columns: k, a.1, a.2 q.1, q.2, level=level, where k is the optimal model, (a.1,a.2) is the interval for alpha in the GIC, (q.1, q.2) is the interval for q and level is the probability. Each row corresponds to an entry in 'level'. |
table |
This table indicates which models can be selected for some values of alpha or q. |
AIC corresponds to setting level=0.84. BIC corresponds to setting level=pchisq(log(n), 1). So for n=100, 1000; BIC=0.96, 0.97
Changjiang Xu and A. Ian McLeod
Changjiang Xu and A. Ian McLeod (2010). Bayesian information criterion with Bernoulli prior. Submitted for publication.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 | #Example 1.
#AR(p) Order selection for 'lynx' series
z <- log(lynx)
n <- length(z)
lag.max <- 20
zta<-ARToPacf(ar.burg(z,aic=FALSE,order.max=lag.max)$ar)
LagsEntering<-1:lag.max
LLapprox <- (-n)*log(cumprod(1-zta[LagsEntering]^2))
ans<-BICqLL(logL=LLapprox, n=n, level=c(0.9, 0.95, 0.99))
ans$khat
ans$table
#if we just want the best model for level=0.99 then,
(BICqLL(logL=LLapprox, n=n, level=0.99)$khat)[[1]]
#aic for comparison
aic<-(-2*LLapprox)+2*LagsEntering
which.min(aic)
plot(LagsEntering, aic)
#
#Example 2. AR(p) Order Selection
#white noise. We do NumRep simulations and
# count the number of overfit models.
set.seed(231789) #make reproducible
n <- 100
lag.max <- 30
LagsEntering<-0:lag.max
NumRep<-25
level<-c(0.99, 0.95, 0.9)
k<-numeric(length(level))
for (i in 1:NumRep){
z <- rnorm(n)
zta<-ARToPacf(ar.burg(z,aic=FALSE,order.max=lag.max)$ar)
LLapprox <- c(0,(-n)*log(cumprod(1-zta[LagsEntering]^2)))
k<-k+as.numeric(0<(BICqLL(logL=LLapprox, n=n, level=level,mSize=LagsEntering)$khat)[,1])
}
ans<-k
names(ans)<-level
ans
#
#Example 3. AR(p) best subset. ARz Family.
z <- log(lynx)
n <- length(z)
lag.max <- 20
zta <- ARToPacf(ar.burg(z,aic=FALSE,order.max=lag.max)$ar)
LagsEntering <- order(abs(zta),decreasing=TRUE)
LLapprox <- c(0, -n*log(cumprod(1-zta[LagsEntering]^2)))
kHat <- (BICqLL(logL=LLapprox, n=n, level=0.99)$khat)[[1]]
pvec<-LagsEntering[1:kHat]
pvec
#pvec above shows the lags in order of importance
#
#Example 4. AR(p) best subset. ARp Family.
#could also try z <- sunspot.year
z <- log(lynx)
lag.max <- 15
pvec <- 1:lag.max
n <- length(z)-lag.max
ind <- (lag.max+1):length(z)
y<-z[ind]
X<-matrix(rep(0,n*lag.max), nrow=n, ncol=lag.max)
for (i in 1:lag.max)
X[,i] <- z[ind-pvec[i]]
outLeaps <- leaps(y=y,x=X,nbest=1,method="r2",strictly.compatible=FALSE)
# approximate likelihood approach
TotSS <- sum((y-mean(y))^2)
RSS <- TotSS*(1-outLeaps$r2)
LogL <- (-n/2)*log(c(TotSS/n, RSS/n))#null model included
ans<-BICqLL(logL=LogL, n=n, level=0.99)
kHat <- (ans$khat)[[1]]-1 #kHat=0 is null model
pvec <-0
if (kHat > 0)
pvec <- (1:lag.max)[(outLeaps$which)[kHat,]]
pvec
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.