r软件mcmc方法

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

我在尝试做一个用mcmc方法寻找变点文章用R做模拟，但我的程序老是运行不对，请各位大神指点指点（本人R菜鸟）结果老是显示这样：Error in if (runif(1) < alpha1) { : missing value where TRUE/FALSE neededIn addition: Warning message:
In rbeta(1, s2 + 1, m * (ck2 - ck1) - s2 + 1) : NAs produced
下面是我的程序

set.seed(2)
x1<-rbinom(8,30,0.3)
x2<-rbinom(6,30,0.8)
x3<-rbinom(6,30,0.6)
y1<-rbinom(8,30,0.75)
y2<-rbinom(6,30,0.75)
y3<-rbinom(6,30,0.75)
c<-rbinom(8,1,0.2)
d<-which(c==1)
x1[d]
T=vector()
for(i in 1:length(d))
T<-qbinom(pbinom(0,30,0.3)+runif(1,min=0,max=pbinom(x1[d],30,0.3)-pbinom(0,30,0.3)),30,0.3)
T
x1[d]<-T
x1
####x1的值

x2
z1=vector()
for(i in 1:6){
z1=ifelse(x2>y2,x2,0)}
f<-which(z1!=0)
T2<-x2[f]
x=vector()
for(i in 1:length(f)){
x=qbinom(pbinom(T2,30,0.3)+runif(1,min=0,max=pbinom(30,30,0.3)-pbinom(T2,30,0.3)),30,0.3)}
x2[f]<-x
f2<-which(z1==0)
T3<-x2[f2]
c2<-rbinom(length(f2),1,0.2)
d2<-which(c2==1)
h=vector()
for(i in 1:length(d2)){
h=qbinom(pbinom(0,30,0.3)+runif(1,min=0,max=pbinom(T3[d2],30,0.3)-pbinom(0,30,0.3)),30,0.3)}
x2[f2][d2]<-h
x2
######x2的值

z3=vector()
for(i in 1:6){
z3=ifelse(x3>y3,x3,0)}
f3<-which(z3!=0)
f4<-which(z3==0)
c3<-rbinom(length(f4),1,0.2)
c3
d3<-which(c3==1)
h2=vector()
for(i in 1:length(d3)){
h2=qbinom(pbinom(0,30,0.3)+runif(1,min=0,max=pbinom(x3[f4][d3],30,0.3)-pbinom(0,30,0.3)),30,0.3)}
x3=x3
x<-c(x1,x2,x3)
s1=sum(x1)
s2=sum(x2)
s3=sum(x3)
m=30
n=20
mchain=matrix(NA,5,1000)
mchain[,1]=c(0.2,0.4,0.6,5,10)
for(i in 2:1000)
{
cp1=mchain[1,i-1]
cp2=mchain[2,i-1]
cp3=mchain[3,i-1]
ck1=mchain[4,i-1]
ck2=mchain[5,i-1]
cp1=rbeta(1,s1+1,m*ck1-s1+1)
cp2=rbeta(1,s2+1,m*(ck2-ck1)-s2+1)
cp3=rbeta(1,s3+1,m*(n-ck2)-s3+1)
pk1=sample(x=seq(2,n-1),size=1)
MHratio1=((1-cp1)*(1-cp2)^(-1))^(m*(pk1-ck1))
alpha1=min(1,MHratio1)
if(runif(1)<alpha1)
{
ck1=pk1
ck1=ck1
}
pk2=sample(x=seq(pk1,n-1),size=1)
MHratio2=((1-cp2)*(1-cp3)^(-1))^(m*(pk2-ck2))
alpha2=min(1,MHratio2)
if(runif(1)<alpha2)
{
ck2=pk2
ck2=ck2
}
mchain[,i]=c(cp1,cp2,cp3,ck1,ck2)
}

扫码加我 拉你入群

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

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

分享0 收藏2 回帖

关键词：mcmc CMC r软件 produced addition

谢谢分享