我要做用mcmc方法求正态分布变点，结果程序老是运行不了，有没有人能给看看，下面是程序：
y1<-rnorm(8,0.8,2)
y2<-rnorm(6,1.2,2)
y3<-rnorm(6,0.4,2)
y<-c(y1,y2,y3)
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<-function(a,b,c)
{
sum(((y-c)^2)[a:b])
}
for(i in 2:numit)
{
y1hat<-sum(y[1:ck1])/ck1
y2hat<-sum(y[(ck1+1):ck2])/(ck2-ck1)
y3hat<-sum(y[(ck2+1):n])/(n-ck2)
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<-exp((like1(1,pk1,ctheta1)+like1(pk1+1,ck2,ctheta2))/(-2*csigma))/exp((like1(1,ck1,ctheta1)+like1(ck1+1,ck2,ctheta2))/(-2*csigma))
(alpha1=min(1,MHratio1))
(ck1<-ifelse((runif(1))<alpha1,pk1,ck1))
(pk2=sample(x=seq(ck1+1,n-1),1))
(MHratio2<-exp((like1(ck1+1,pk2,ctheta2)+like1(pk2+1,n,ctheta3))/(-2*csigma))/exp((like1(ck1+1,ck2,ctheta2)+like1(ck2+1,n,ctheta3))/(-2*csigma))
(alpha2=min(1,MHratio2))
(ck2<-ifelse((runif(1))<alpha2,pk2,ck2))
mchain[,i]=c(ctheta1,ctheta2,ctheta3,csigma,ck1,ck2)
}
return(mchain)
}
apply(mhsampler(1000,dat=y),1,mean)