全部版块 我的主页
论坛 提问 悬赏 求职 新闻 读书 功能一区 经管百科 爱问频道
873 1
2018-12-03
j<-0;u0<-5;u10<-12;u11<-10;o0<-4;o10<-5;o11<-8;p<-0.1
repeat{
  j<-j+1
  q1<-0;q2<-0;q3<-0;q4<-0;q5<-0;q6<-0;q7<-0;q8<-0;q9<-0
  for(i in 1:m1){
    f1<-function(a,b){ #qab
      v1<-(1/k)*(a*u10+b*u11+(k-a-b)*u0)#uab'
      c1<-(1/k^2)*(a*o10+b*o11+(k-a-b)*o0)#oab'
      s1<-dbinom(a,e[i],p)*dbinom(b,k-e[i],p)*(2*pi*c1)^(-1/2)*exp(-(d[i]-v1)^2/(2*c1))
      s2<-0
      for(a2 in 0:e[i]){
        for(b2 in 0:k-e[i]){
          v2<-(1/k)*(a2*u10+b2*u11+(k-a2-b2)*u0)#uab'
          c2<-(1/k^2)*(a2*o10+b2*o11+(k-a2-b2)*o0)#oab'
          s2<-s2+dbinom(a2,e[i],p)*dbinom(b2,k-e[i],p)*(2*pi*c2)^(-1/2)*exp(-(d[i]-v2)^2/(2*c2))
        }
      }
      return(s1/s2)
    }
    f2<-function(a,b){ #qab*(yi-uab')^2
      v1<-(1/k)*(a*u10+b*u11+(k-a-b)*u0)#uab'
      c1<-(1/k^2)*(a*o10+b*o11+(k-a-b)*o0)#oab'
      s3<-dbinom(a,e[i],p)*dbinom(b,k-e[i],p)*(2*pi*c1)^(-1/2)*exp(-(d[i]-v1)^2/(2*c1))*(d[i]-v1)^2
      s4<-0
      for(a2 in 0:e[i]){
        for(b2 in 0:k-e[i]){
          v2<-(1/k)*(a2*u10+b2*u11+(k-a2-b2)*u0)#uab'
          c2<-(1/k^2)*(a2*o10+b2*o11+(k-a2-b2)*o0)#oab'
          s4<-s4+dbinom(a2,e[i],p)*dbinom(b2,k-e[i],p)*(2*pi*c2)^(-1/2)*exp(-(d[i]-v2)^2/(2*c2))
        }
      }
      return(s3/s4)
    }
    q1<-q1+f1(0,0)*d[i]#q00*yi连加
    q2<-q2+f1(0,0)#q00连加
    q3<-q3+f1(0,1)*d[i]#q01*yi连加
    q4<-q4+f1(0,1)#q01连加
    q5<-q5+f1(1,0)*d[i]#q10*yi连加
    q6<-q6+f1(1,0)#q10连加
    q7<-q7+f2(0,0)#q00*(yi-u00')^2连加
    q8<-q8+f2(0,1)#q01*(yi-u01')^2连加
    q9<-q9+f2(1,0)#q10*(yi-u10')^2连加
  }
  oldu0<-u0;oldu11<-u11;oldu10<-u10;oldo0<-o0;oldo11<-o11;oldo10<-o10;oldp<-p#记录上一次的迭代值
  u0<-q1/q2
  u11<-k*q3/q4-(k-1)*u0
  u10<-k*q5/q6-(k-1)*u0
  o0<-k*q7/q2
  o11<-(k^2)*q8/q4-(k-1)*o0
  o10<-(k^2)*q9/q6-(k-1)*o0
  p<-1-(q2/m1)^(1/k)
  epsilo<-1
  if((u0-oldu0)^2<epsilo&
     (u11-oldu11)^2<epsilo&
     (u10-oldu10)^2<epsilo&
     (o0-oldo0)^2<epsilo&
     (o10-oldo10)^2<epsilo&
     (o11-oldo11)^2<epsilo&
     (p-oldp)^2<epsilo)break
}
Error in if ((u0 - oldu0)^2 < epsilo & (u11 - oldu11)^2 < epsilo & (u10 -  :
  missing value where TRUE/FALSE needed
> u0
[1] NaN
> u11
[1] NaN

二维码

扫码加我 拉你入群

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

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

全部回复
2018-12-4 14:22:33
您好,如果您的求助没有解决,请到项目交易发布需求,会有更快更专业的用户帮助您 https://bbs.pinggu.org/prj/

项目交易是为用户提供需求的平台,可以在平台发布你需求,也可以展现你的技术帮助他人,从而得到相应的报酬。
二维码

扫码加我 拉你入群

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

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

栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群