Мои вопросы касаются сравнения двух методов прогнозирования биномиального исхода после подгонки модели bsts к данным обучения.Метод A использует предикат.bsts, а B напрямую использует конечное состояние и подогнанные коэффициенты регрессии.B, по-видимому, дает более точные прогнозы и дает выборки из прогнозируемой вероятности (включающие CI), в отличие от A, который обеспечивает биномиальную ничью.Вопросы: неправильно ли я понял A, почему методы дают разные результаты, и действительно ли B лучше?
Ниже A и B используются для прогнозирования ответа в момент времени j на основе ответов и ковариат в моменты времени 1:(j-1) вместе с ковариатами в момент времени j.Априорные значения в соответствии с примером 3 в http://www.unofficialgoogledatascience.com/2017/07/fitting-bayesian-structural-time-series.html
set.seed(1)
x1<-rnorm(100)
x2<-rpois(100,3)
x3<-runif(100)
g<-function(t) sin(pi*t/100)
lp<-2*g(1:100)+1.8*x1-x2
#trend, and no dependence on x3
h<-function(t) (1+exp(-t))^-1
p<-h(lp)
y<-rbinom(N,1,p)
burn=1:1000
ss<-AddLocalLevel(list(),
sigma.prior=SdPrior(sigma.guess=0.1,sample.size=1,upper.limit=1),
initial.state.prior=NormalPrior(0,5))
mod1<-bsts(y~-1+x1+x2+x3,family="logit",
state.specification=ss,niter=2000,seed=12345)
plot(mod1,"comp")
plot(mod1,"coef")
S=10
fit1<-fit2<-rep(NA,N)
pb<-txtProgressBar(min=S,max=N,initial=S,style=3)
for (j in S:N)
{
setTxtProgressBar(pb,value=j)
yj<-y[1:(j-1)]
x1j<-x1[1:(j-1)]
x2j<-x2[1:(j-1)]
x3j<-x3[1:(j-1)]
ssj<-AddLocalLevel(list(),
sigma.prior=SdPrior(sigma.guess=0.1,
sample.size=1,upper.limit=1),
initial.state.prior=NormalPrior(0,5))
mod1j<-bsts(yj~-1+x1j+x2j+x3j,
family="logit",state.specification=ssj,
niter=2000,seed=12345)
newj<-data.frame(x1j=x1[j],x2j=x2[j],x3j=x3[j])
#method A
predj<-predict(mod1j,newdata=newj,seed=12345)
fit1[j]<-predj$mean
#method B
preddj<-sapply(mod1j$fin[-burn,]+mod1j$coef[-burn,]%*%t(as.matrix(newj)),h)
fit2[j]<-mean(preddj)
}
mse1<-round(sum((y[S:N]-fit1[S:N])^2)/(N-S+1),2)
mse2<-round(sum((y[S:N]-fit2[S:N])^2)/(N-S+1),2)
mse1x<-round(sum((y[(N/2+1):N]-fit1[(N/2+1):N])^2)/(N/2),2)
mse2x<-round(sum((y[(N/2+1):N]-fit2[(N/2+1):N])^2)/(N/2),2)
plot(1:N,y,col=ifelse(y==1,"red","blue"))
for (j in 1:N)
{
points(j,fit1[j],pch=20,col="seagreen")
points(j,fit2[j],col="brown")
lines(c(j,j),c(0,1),lty=2,col="grey")
}
text(0,0.5,pos=4,paste("mse1=",mse1,sep=""))
text(0,0.45,pos=4,paste("mse2=",mse2,sep=""))
text(0,0.3,pos=4,paste("mse1x=",mse1x,sep=""))
text(0,0.25,pos=4,paste("mse2x=",mse2x,sep=""))
#mse2<mse1 and mse2x<mse1x
моделирования
В 50 наборах смоделированных данных коэффициенты, примененные к g, x1 и x2, были установлены с помощью rnorm (3).В 41 прогоне mse2 ниже mse1.В 42 прогонах mse2x ниже mse1x.Средние значения mse1, mse2, mse1x и mse2x составили 0,2366, 0,1418, 0,2426, 0,1354 соответственно.