Сообщение об ошибке в алгоритме EM для оценки параметра масштаба - PullRequest
0 голосов
/ 11 октября 2018

Это мой код алгоритма EM для оценки масштабного параметра кругового распределения.Сначала я имитирую данные из двумерной нормали, а затем преобразую в углы, чтобы получить круговое распределение.Параметр масштаба, записанный в функции: (varrr).Соответствующая оценка (EMest2).

library(OpenMx) ## to use vec2diag command
library (MASS)  ## Simulate from a Multivariate Normal Distribution
n=10 # sample size
eps=10^-8#stopping value
mu=c(1,0.9425) # true values of the mean vector of bivariate normal 
variables
var.1 =1 # sigma.1.squared for first normal variable
var.2=1#sigma.2.squared of second normal variable
cov.1.2=0 # covariance between the two normal variables
var.cov<-matrix(c(var.1,cov.1.2,cov.1.2,var.2),nrow=2,ncol=2,byrow=T)#variance covariance matrix of the bivariate normal random variables    
y= mvrnorm(n, mu, var.cov)#generating the bivariate normal random variables
theta.rad=rep(NA,times=n)
  for (ii in 1:n){
    if (y[ii,1] > 0 & y[ii,2]>=0) {theta.rad[ii]=atan(y[ii,2]/y[ii,1])}
    else if (y[ii,1] < 0) {theta.rad[ii]=(atan(y[ii,2]/y[ii,1]))+pi}
    else if (y[ii,1] > 0 & y[ii,2]<0) {theta.rad[ii]=atan(y[ii,2]/y[ii,1])+2*pi}
    else if (y[ii,1] == 0 & y[ii,2]>0) {theta.rad[ii]=pi/2}
    else if (y[ii,1] ==0 & y[ii,2]<0) {theta.rad[ii]=3*pi/2}
  }# calculating the angles
   u=cbind(cos(theta.rad),sin(theta.rad)) #constructing the main matrix
  var.1.initial=var(y[,1])
  EM.est=function(uu,varrr){
    tta=as.vector(1/sqrt(varrr))
    tt=tta*uu%*%mu
    phidash=(pnorm(tt)/((tt*pnorm(tt))+dnorm(tt)))+tt  
    m=vec2diag(phidash)
    aa=2*varrr-((sum(m%*%tt))/n)+(t(mu)%*%mu)
    return(aa)
  }
  w=0      ###counter 
  diff=1
  while(diff>eps){ 
    EMest1=var.1.initial
    EMest2=EM.est(u,EMest1)
    var.1.initial=EMest2
    diff=max(abs((EMest2-EMest1)/EMest1))
    w=w+1
  }
  print(EMest2)

Каждый раз, когда я запускаю программу, я получаю сообщение об ошибке: Error in while (diff > eps) {:missing value where TRUE/FALSE needed и значение оценки (EMest2) равно бесконечности.Собственно, я не знаю, в чем заключается ошибка в написанной программе.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...