mchain r语言_R语言mcmc方法

我要做用mcmc方法求正态分布变点，结果程序老是运行不了，有没有人能给看看，下面是程序：

y1

y2

y3

y

n=20

ck1=5

ck2=10

mhsampler=function(numit,dat=y)

{

mchain=matrix(NA,6,numit)

mchain[,1]=c(0.2,0.4,0.6,1,5,10)

like1

{

sum(((y-c)^2)[a:b])

}

for(i in 2:numit)

{

y1hat

y2hat

y3hat

ctheta1=mchain[1,i-1]

ctheta2=mchain[2,i-1]

ctheta3=mchain[3,i-1]

csigma=mchain[4,i-1]

ck1=mchain[5,i-1]

ck2=mchain[6,i-1]

(ctheta1=rnorm(1,y1hat,csigma/ck1))

(ctheta2=rnorm(1,y2hat,csigma/(ck2-ck1)))

(ctheta3=rnorm(1,y3hat,csigma/(n-ck2)))

csigma=1/rgamma(1,n-1,(like1(1,ck1,y1hat)+like1(ck1+1,ck2,y2hat)+like1(ck2+1,n,y3hat))/2)

(pk1=sample(x=seq(2,ck2-1),1))

(MHratio1

(alpha1=min(1,MHratio1))

(ck1

(pk2=sample(x=seq(ck1+1,n-1),1))

(MHratio2

(alpha2=min(1,MHratio2))

(ck2

mchain[,i]=c(ctheta1,ctheta2,ctheta3,csigma,ck1,ck2)

}

return(mchain)

}

apply(mhsampler(1000,dat=y),1,mean)