第一次发帖。。能不能帮我看下广义极值分布这个MCMC哪里有错误

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

current.cita<-c(0.1,0.1,0.1)
> cita<-matrix(nrow=m,ncol=3)  
> m<-55000                     
> mu.cita<-c(0,0,0)
> s.cita<-c(10000,10000,100)
> prop.s<-c(0.1,0.1,0.1)
> cita<-matrix(nrow=m,ncol=2)
> acc.prob<-0
> for(t in 1:m){
+ prop.cita<-rnorm(3,current.cita,prop.s)
+  h<-1-current.cita[3]/exp(current.cita[2])*(x-current.cita[1])
g<-1-prop.cita[3]/exp(prop.cita[2])*(x-prop.cita[1])
+ if(prod(h>0)==0){loga<-0}
+ else{loga<--770*current.cita[2]+(1/current.cita[3]-1)*sum(log(1-current.cita[3]/exp(current.cita[2])*(x-current.cita[1])))
+ }
+ if(prod(g>0)==0){logb<-0}
+ else{logb<--770*prop.cita[2]+(1/prop.cita[3]-1)*sum(log(1-prop.cita[3]/exp(prop.cita[2])*(x-prop.cita[1])))
+ }
+ if(prod(g>0)*prod(h>0)==0){logc<-10}
+logc<-sum(dnorm(prop.cita,mu.cita,sigma.cita,log=T))-sum(dnorm(current.cita,mu.cita,sigma.cita,log=T))
+ logl<-logb-loga+logc
+ u<-runif(1)
+ u<-log(u)
+ if(u<logl){
+ current.cita<-prop.cita
+ acc.prob<-acc.prob+1
+ }
+ cita[t,]<-current.cita
+ }

扫码加我 拉你入群

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

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

分享0 收藏1 回帖

关键词：第一次发帖 mcmc 极值分布 CMC 哪里有

你上面定义的是s.cita，下面出现的是sigma.cita
另外，cita[t,]<-current.cita 两个长度都不一样
处理成cita[]就不报错了