Некоторые основные вопросы - не о том, как запустить программу, а о том, что она на самом деле вычисляет. Я не могу получить ответы, набрав bsts, что приводит к Call ("analysis_common_r_fit_bsts_model _" ...), к которому я не могу получить доступ и, возможно, не пойму, если смогу. Я попробовал r-help и отправил электронное письмо сопровождающему. Я надеюсь, что один из вас может помочь или указать мне на ответ.
Для случая Гаусса без регрессоров и модели локального уровня приведенный ниже код, по-видимому, имитирует bsts. Он создает графики, аналогичные графику модели, и подтверждает связь между выходными данными sigma.level и sigma.obs с вкладом состояния. Тем не менее, прямой подход к заднему распределению для этих сигм не совпадает с выходом bsts. Я не уверен, что это простая ошибка, если я неправильно понимаю SdPrior или различия связаны с mcmc ...
setwd("/home/greg/Documents/rwork/incid/bayes")
library(bsts)
library(mcmc)
#simulate some data
y<-rep(NA,50)
y[1]=1
y[2]=1
s=2
set.seed(5)
for (k in 1:48)
{
y[2+k]=y[1+k]+0.1*y[k]+s*rnorm(1)
}
plot(1:50,y[1:50],main=paste("seed =",5))
#bsts model
ss<-AddLocalLevel(list(),y)
mod1<-bsts(y,state.specification=ss,niter=1000)
plot(mod1)
#what is actually being plotted?
#reproduce the plot using mod1$state.contributions
par(mfrow=c(1,2))
plot(1:50,y,col="blue",ylim=c(-12,8),main="quantiles of mod1$state.contributions")
for (i in 1:99)
{
qi<-qin<-rep(NA,50)
tj<-1:50
for (j in 1:50)
{
qi[j]<-quantile(mod1$state.con[501:1000,1,j],(i-0.5)/100)
qin[j]<-quantile(mod1$state.con[501:1000,1,j],(i+1-0.5)/100)
}
polygon(c(tj,rev(tj)),c(qi[1:50],rev(qin[1:50])),col=rgb(0,0,0,40*dnorm(i,mean=50,sd=20)),border=FALSE)
}
plot(mod1,ylim=c(-12,8),main="plot(mod1)")
par(mfrow=c(1,1))
#the state specification is based on sd(y)
sd(y)
ss
#kalman iteration, then add P error term when generating state from xe
#set P[1,,] = 1 to start
#then at end, reset P[1,,]=P[2,,]
#polygon shading
ax<-(mod1$sigma.level[501:1000])^2
bx<-(mod1$sigma.obs[501:1000])^2
state<-matrix(NA,500,50)
for (j in 1:500)
{
a<-ax[j]
b<-bx[j]
H=matrix(1,1,1)
F=matrix(1,1,1)
N=50
dim(y)=c(1,N)
xe<-ye<-ze<-matrix(NA,1,N)
xe[,1]<-1
ye[,1]<-H%*%xe[,1]
ze[,1]<-y[,1]-ye[,1]
P<-K<-array(data=NA,dim=c(N,1,1))
P[1,,]<-1
for (i in 1:(N-1))
{
P[i+1,,]<-F%*%P[i,,]%*%t(F)+a
K[i+1,,]<-P[i+1,,]%*%t(H)%*%solve(b+H%*%P[i+1,,]%*%t(H))
P[i+1,,]<-(diag(1,1)-K[i+1,,]%*%H)%*%P[i+1,,]
xe[1,i+1]<-F%*%xe[,i]+K[i+1,,]%*%(y[,i+1]-H%*%F%*%xe[,i])
}
P[1,,]<-P[2,,]
state[j,]<-rnorm(N,xe[1,1:N],P[1:N,,]^0.5)
}
par(mfrow=c(1,2))
plot(1:N,y,col="blue",ylim=c(-12,8),main="quantiles of kalman state")
for (i in 1:99)
{
qi<-qin<-rep(NA,N)
tj<-1:N
for (k in 1:N)
{
qi[k]<-quantile(state[,k],(i-0.5)/100)
qin[k]<-quantile(state[,k],(i+1-0.5)/100)
}
polygon(c(tj,rev(tj)),c(qi[1:N],rev(qin[1:N])),col=rgb(0,0,0,40*dnorm(i,mean=50,sd=20)),border=FALSE)
}
plot(mod1,ylim=c(-12,8))
par(mfrow=c(1,1))
#indices chosen at random
#plausible that state[,i] and mod1$state.con[501:1000,1,i] are the same distribution
par(mfrow=c(2,4),mar=c(5,5,3,1))
sami<-sample(1:50,8,replace=F)
for (i in 1:8)
{
qqplot(state[,sami[i]],mod1$state.con[501:1000,1,sami[i]],cex=0.5,main=paste("i =",sami[i]),cex.lab=1.3,xlab=paste("state[,",sami[i],"]",sep=""),ylab=paste("mod1$state.con[501:1000,1,",sami[i],"]",sep=""))
lines(state[,sami[i]],state[,sami[i]],col="red")
}
par(mfrow=c(1,1))
#the likelihood, using the kalman filter, as a function of the error variances and the initial state
kal<-function(par)
{
a<-par[1]
b<-par[2]
init<-par[3]
initP<-par[4]
H=matrix(1,1,1)
F=matrix(1,1,1)
#1-dimensional state
N=50
dim(y)=c(1,N)
xe<-ye<-matrix(NA,1,N)
xe[,1]<-init
ye[,1]<-H%*%xe[,1]
P<-K<-array(data=NA,dim=c(N,1,1))
#P[1,,] initial guess
P[1,,]<-initP
for (i in 1:(N-1))
{
P[i+1,,]<-F%*%P[i,,]%*%t(F)+a
K[i+1,,]<-P[i+1,,]%*%t(H)%*%solve(b+H%*%P[i+1,,]%*%t(H))
xe[1,i+1]<-F%*%xe[,i]+K[i+1,,]%*%(y[,i+1]-H%*%F%*%xe[,i])
P[i+1,,]<-(diag(1,1)-K[i+1,,]%*%H)%*%P[i+1,,]
}
-1/2*(log(abs(b))+(1/b)*sum((y[1,]-xe[1,])^2))
}
#independent priors on a and b
#bsts uses sd(y) to set both priors.
#for a (ss[[1]]$sigma.prior) it uses SdPrior with
#$prior.guess 0.04600655, $prior.df 0.01, $initial.value 0.04600655, $upper.limit 4.600655
#for b (mod1$prior) it uses SdPrior with
#$prior.guess 4.600655, $prior.df 0.01, $initial.value 4.600655, $upper.limit 5.520786
#ss[[1]]$initial.state.prior is normal with mu=1, sd=4.600655
#try an inverse gamma prior on a and b, normal prior on init, with initP fixed
initP<-0.25
lpr1<-function(a)
{
v=0.01
ifelse((a<=0)|a>(sd(y))^2,-Inf,-(v/2+1)*log(a)-v*(sd(y)/100)^2/(2*a))
}
lpr2<-function(b)
{
v=0.01
ifelse((b<=0)|b>(1.2*sd(y))^2,-Inf,-(v/2+1)*log(b)-v*(sd(y))^2/(2*b))
}
lpr3<-function(c)
{
-1/(2*sd(y))*(c-1)^2
}
lpost1<-function(par)
{
a<-par[1]
b<-par[2]
init<-par[3]
lpr1(a)+lpr2(b)+lpr3(init)+kal(c(par,initP))
}
par0<-c(sd(y)/100,sd(y),1)
nb=5000
out<-metrop(lpost1,par0,nb,scale=0.5)
sam<-out$batch
dim(sam)
print(sam[500:510,],12)
par(mfrow=c(3,1))
plot(sam[(nb/2):nb,1],type="l")
plot(sam[(nb/2):nb,2],type="l")
plot(sam[(nb/2):nb,3],type="l")
par(mfrow=c(1,1))
op<-optim(par0,lpost1,control=list(fnscale=-1))
out<-metrop(lpost1,op$par,nb,scale=1)
sam<-out$batch
dim(sam)
print(sam[500:510,],12)
par(mfrow=c(3,1))
plot(sam[(nb/2):nb,1],type="l")
plot(sam[(nb/2):nb,2],type="l")
plot(sam[(nb/2):nb,3],type="l")
par(mfrow=c(1,1))
mpar<-apply(sam[(nb/2):nb,],2,mean)
lpost1m<-function(par) lpost1(par+mpar)
opm<-optim(par0,lpost1m,control=list(fnscale=-1))
#a little fiddling gets plausible mixing
nb=5000
#be patient
outm<-metrop(lpost1m,opm$par,nb,blen=5,nspac=10,scale=1.5)
outm$acc
samm<-outm$batch
dim(samm)
par(mfrow=c(3,1))
plot(samm[(nb/2):nb,1],type="l")
plot(samm[(nb/2):nb,2],type="l")
plot(samm[(nb/2):nb,3],type="l")
par(mfrow=c(1,1))
acf(samm)
samm<-thin(samm,5)
dim(samm)
acf(samm)
print(samm[501:510,],8)
par(mfrow=c(3,1))
plot(samm[501:1000,1],type="l")
plot(samm[501:1000,2],type="l")
plot(samm[501:1000,3],type="l")
par(mfrow=c(1,1))
#plot kalman using samm from lpost1m
ax<-samm[501:1000,1]+mpar[1]
bx<-samm[501:1000,2]+mpar[2]
initx<-samm[501:1000,3]+mpar[3]
state<-matrix(NA,500,50)
for (j in 1:500)
{
a<-ax[j]
b<-bx[j]
init<-initx[j]
H=matrix(1,1,1)
F=matrix(1,1,1)
N=50
dim(y)=c(1,N)
xe<-ye<-matrix(NA,1,N)
xe[,1]<-init
ye[,1]<-H%*%xe[,1]
P<-K<-array(data=NA,dim=c(N,1,1))
P[1,,]<-0
for (i in 1:(N-1))
{
P[i+1,,]<-F%*%P[i,,]%*%t(F)+a
K[i+1,,]<-P[i+1,,]%*%t(H)%*%solve(b+H%*%P[i+1,,]%*%t(H))
P[i+1,,]<-(diag(1,1)-K[i+1,,]%*%H)%*%P[i+1,,]
xe[1,i+1]<-F%*%xe[,i]+K[i+1,,]%*%(y[,i+1]-H%*%F%*%xe[,i])
}
P[1,,]<-P[2,,]
state[j,]<-rnorm(N,xe[1,1:N],P[1:N,,]^0.5)
}
par(mfrow=c(1,2))
plot(1:N,y,col="blue",ylim=c(-12,8))
for (i in 1:99)
{
qi<-qin<-rep(NA,N)
tj<-1:N
for (k in 1:N)
{
qi[k]<-quantile(state[,k],(i-0.5)/100)
qin[k]<-quantile(state[,k],(i+1-0.5)/100)
}
polygon(c(tj,rev(tj)),c(qi[1:N],rev(qin[1:N])),col=rgb(0,0,0,35*dnorm(i,mean=50,sd=15)),border=FALSE)
}
plot(mod1,ylim=c(-12,8))
par(mfrow=c(1,1))
#plausible
#but
mean(ax)
mean((mod1$sigma.level[501:1000])^2)
qqplot(ax,(mod1$sigma.level[501:1000])^2)
mean(bx)
mean((mod1$sigma.obs[501:1000])^2)
qqplot(bx,(mod1$sigma.obs[501:1000])^2)
lines(bx,bx,col="red")
#so ax is not sampling (sigma.level)^2
#bx not sampling (sigma.obs)^2