удаляя итерации, которые не удовлетворяли условию в R

Question

Я выполняю симуляцию для 10000 прогонов модели местоположения.проблема в том, что когда я вычисляю оценку m (EM.var.est), я хотел бы удалить итерации с большими значениями (> 100) и заменить их другими итерациями.счетчик для номера итерации равен (i), я ставлю (i = i-1) в конце цикла, так что если условие не выполняется (i) не учитывается, но кажется, что программа этого не делает, иЯ получил результат NA в окончательном векторе (vec.var).Вот пример кода

    #---- Prepare the R console
rm(list = ls(all = TRUE)) # Remove all objects in R console
set.seed(123456) # Set the seed for reproducible results.default for the program to fix initial values to start with in each time.
#----------------------
library (MASS)  ## Simulate from a Multivariate Normal Distribution
runs=10000 #number of repeated samples
n=10 # sample size
mu=c(0.1,0.3816) # true values of the mean vector of bivariate normal variables
var.1 =0.1# sigma.1.squared for first normal variable
var.2=0.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    
vec.var=rep(NA,times=runs)
for (i in 1:runs){
  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=apply(u,2,var)[1] # initial value for the variance
  EM.est=function(uu,varrr){
    tta=1/sqrt(varrr)
    tt=tta*uu%*%mu
    phidash=(pnorm(tt)/((tt*pnorm(tt))+dnorm(tt)))+tt  
    m=diag(drop(phidash))
    aa=drop(varrr*(2-((sum(m%*%tt))/n))+(t(mu)%*%mu))
    return(aa)
  }
  w=0 ###counter 
  ll1=0
  ll2=1
  while(ll1<ll2){ 
    EMest1=var.1.initial
    a1=(-n*(log(2*pi)))
    a2=((-n/(2*EMest1))*(t(mu)%*%mu))
    q=(1/sqrt(EMest1))*u%*%mu
    a3=log(1+((q*pnorm(q))/dnorm(q)))
    ll1=a1+a2+sum(a3)
    EMest2=EM.est(u,EMest1)
    a21=(-n*(log(2*pi)))
    a22=((-n/(2*EMest2))*(t(mu)%*%mu))
    q2=(1/sqrt(EMest2))*u%*%mu
    a23=log(1+((q2*pnorm(q2))/dnorm(q2)))
    ll2=a21+a22+sum(a23)
    var.1.initial=EMest2
    w=w+1
  }
  EM.est.var=EMest2
  if (EM.est.var<100) {vec.var[i]=EM.est.var}
  else {i=i-1}
}

summary(vec.var)

Мой вопрос заключается в том, как я могу сделать еще одну итерацию вместо этой и сохранить итоговое количество итераций, как есть. Большое спасибо

Answer 1

Вы можете использовать цикл while вместо for и увеличивать значение i только при соблюдении условия.Это гарантирует, что симуляция выдает runs действительных значений.

rm(list = ls(all = TRUE)) # Remove all objects in R console
set.seed(123456) # Set the seed for reproducible results.default for the program to fix initial values to start with in each time.
#----------------------
library (MASS)  ## Simulate from a Multivariate Normal Distribution
runs=10000 #number of repeated samples
n=10 # sample size
mu=c(0.1,0.3816) # true values of the mean vector of bivariate normal variables
var.1 =0.1# sigma.1.squared for first normal variable
var.2=0.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    
vec.var=rep(NA,times=runs)
i <- 1
while(i <= runs){
#for (i in 1:runs){
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=apply(u,2,var)[1] # initial value for the variance
EM.est=function(uu,varrr){
    tta=1/sqrt(varrr)
    tt=tta*uu%*%mu
    phidash=(pnorm(tt)/((tt*pnorm(tt))+dnorm(tt)))+tt  
    m=diag(drop(phidash))
    aa=drop(varrr*(2-((sum(m%*%tt))/n))+(t(mu)%*%mu))
    return(aa)
}
w=0 ###counter 
ll1=0
ll2=1
while(ll1<ll2){ 
    EMest1=var.1.initial
    a1=(-n*(log(2*pi)))
    a2=((-n/(2*EMest1))*(t(mu)%*%mu))
    q=(1/sqrt(EMest1))*u%*%mu
    a3=log(1+((q*pnorm(q))/dnorm(q)))
    ll1=a1+a2+sum(a3)
    EMest2=EM.est(u,EMest1)
    a21=(-n*(log(2*pi)))
    a22=((-n/(2*EMest2))*(t(mu)%*%mu))
    q2=(1/sqrt(EMest2))*u%*%mu
    a23=log(1+((q2*pnorm(q2))/dnorm(q2)))
    ll2=a21+a22+sum(a23)
    var.1.initial=EMest2
    w=w+1
}
EM.est.var=EMest2
if (EM.est.var<100) {vec.var[i]=EM.est.var; i <- i+1}
#  else {i=i-1}
}

summary(vec.var)