R语言mcmc方法

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

我要做用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)

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏4 回帖

关键词：mcmc CMC R语言 正态分布 有没有人

不要沉啊！！！！！！！自己顶一下

1. y1<-rnorm(8,0.8,2)
   
2. y2<-rnorm(6,1.2,2)
   
3. y3<-rnorm(6,0.4,2)
   
4. y<-c(y1,y2,y3)
   
5. n=20
   
6. ck1=5
   
7. ck2=10
   
8. mhsampler=function(numit,dat=y)
   
9. {
   
10. mchain=matrix(NA,6,numit)
    
11. mchain[,1]=c(0.2,0.4,0.6,1,5,10)
    
12. like1<-function(a,b,c)
    
13. {
    
14. sum(((y-c)^2)[a:b])
    
15. }
    
16. for(i in 2:numit)
    
17. {
    
18. y1hat<-sum(y[1:ck1])/ck1
    
19. y2hat<-sum(y[(ck1+1):ck2])/(ck2-ck1)
    
20. y3hat<-sum(y[(ck2+1):n])/(n-ck2)
    
21. ctheta1=mchain[1,i-1]
    
22. ctheta2=mchain[2,i-1]
    
23. ctheta3=mchain[3,i-1]
    
24. csigma=mchain[4,i-1]
    
25. ck1=mchain[5,i-1]
    
26. ck2=mchain[6,i-1]
    
27. (ctheta1=rnorm(1,y1hat,csigma/ck1))
    
28. (ctheta2=rnorm(1,y2hat,csigma/(ck2-ck1)))
    
29. (ctheta3=rnorm(1,y3hat,csigma/(n-ck2)))
    
30. csigma=1/rgamma(1,n-1,(like1(1,ck1,y1hat)+like1(ck1+1,ck2,y2hat)+like1(ck2+1,n,y3hat))/2)
    
31. (pk1=sample(x=seq(2,ck2-1),1))
    
32. (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)))
    
33. (alpha1=min(1,MHratio1))
    
34. (ck1<-ifelse((runif(1))<alpha1,pk1,ck1))
    
35. (pk2=sample(x=seq(ck1+1,n-1),1))
    
36. (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)))
    
37. (alpha2=min(1,MHratio2))
    
38. (ck2<-ifelse((runif(1))<alpha2,pk2,ck2))
    
39. mchain[,i]=c(ctheta1,ctheta2,ctheta3,csigma,ck1,ck2)
    
40. }
    
41. return(mchain)
    
42. }
    
43. apply(mhsampler(1000,dat=y),1,mean)

复制代码

MHratio1，MHratio2最后少加了一个括号，你能给我发一份推到过程吗？1769138434@qq.com，谢谢

zhou1_20 发表于 2017-3-11 10:42
MHratio1，MHratio2最后少加了一个括号，你能给我发一份推到过程吗？1769138434@qq.com，谢谢

谢谢！！我最近没逛论坛，在准备开题，我给你发送了

Metropolis-Hasting算法有R的包的，不必自己写～

foozhencheng 发表于 2017-5-28 13:50
Metropolis-Hasting算法有R的包的，不必自己写～

你好，能讲讲用什么包么

Zjn55571993 发表于 2017-6-1 16:28
你好，能讲讲用什么包么

可以查一下MHadaptive

谢谢！

好东西，谢谢分享