bsts 后验分布答案

【问题标题】：bsts posterior distributionsbsts 后验分布
【发布时间】：2019-01-19 11:00:31
【问题描述】：

提出了一些基本问题——不是如何运行程序，而是它实际在计算什么。我无法通过键入 bsts 来获得答案，这会导致 Call("analysis_common_r_fit_bsts_model_"...) 我无法访问并且可能无法理解。我已经尝试过 r-help 并向包维护者发送电子邮件。我希望你们中的一个可以帮助或指出我的答案。

对于没有回归量和局部级别模型的高斯情况，下面的代码似乎模仿了 bsts。它生成类似于模型图的图形，并确认 sigma.level 和 sigma.obs 的输出与状态贡献之间的关系。然而，对这些 sigma 的后验分布的直接方法与 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

【问题讨论】：

标签： r

【解决方案1】：

现在有所改善，主要是通过更改先验中的最大值。

后验均值更接近 bsts 输出，但分布仍然不同。我不知道 bsts 中使用的是什么 mcmc 实现。

格雷格

#

setwd("/home/greg/Documents/rwork/incid/bayes")
library(bsts)
library(mcmc)
library(coda)

#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=5000)
plot(mod1)

#has mod1 converged?
sigmc<-as.mcmc(cbind(mod1$sigma.level,mod1$sigma.obs))
plot(sigmc)
acf(sigmc)
#highly correlated
heidel.diag(sigmc)
geweke.plot(sigmc)
#looks ok

#skipping the discussion of what is being plotted

#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]

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,,]<-b

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))
}


#priors

lpr1<-function(a)
{
v=0.01
ifelse((a<=0)|a>sd(y),-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),-Inf,-(v/2+1)*log(b)-v*(sd(y))^2/(2*b))
}

lpr3<-function(c)
{
-1/(2*(sd(y))^2)*(c-1)^2
}

lpost1<-function(par)
{
a<-par[1]
b<-par[2]
init<-par[3]
lpr1(a)+lpr2(b)+lpr3(init)+kal(par)
}

par0<-c(sd(y)/100,sd(y),1)

nb=5000
out<-metrop(lpost1,par0,nb,blen=5,nspac=5,scale=0.3)
out$acc
sam<-out$batch
samc<-as.mcmc(sam)
plot(samc)
heidel.diag(samc)
geweke.plot(samc)
acf(samc)

#plot kalman using samc

ax<-samc[2501:5000,1]
bx<-samc[2501:5000,2]
initx<-samc[2501:5000,3]

state<-matrix(NA,2500,50)

for (j in 1:2500)
{
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,,]<-b
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])
}

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,40*dnorm(i,mean=50,sd=20)),border=FALSE)
}
plot(mod1,ylim=c(-12,8))
par(mfrow=c(1,1))

#plausible
#but

mean(ax^0.5)
mean(mod1$sigma.level[2501:5000])
#better
qqplot(ax^0.5,mod1$sigma.level[2501:5000])
lines(ax^0.5,ax^0.5,col="red")

mean(bx^0.5)
mean(mod1$sigma.obs[501:1000])
qqplot(bx^0.5,mod1$sigma.obs[501:1000])
lines(bx^0.5,bx^0.5,col="red")

#so ax^0.5 is still not sampling sigma.level
#bx^0.5 still not sampling sigma.obs

【讨论】：