1

A partir de unos datos que estoy generando (weibull) estoy tratando de hacer un Phase Type fitting usando el algortimo EM en R. Sin embargo, al generar la función de densidad (matrizp) me sale el siguiente error:

"matrizp[j,i]= pmf(data[j],lambds[i],rs[i]) replacement has length
zero"

El código es el siguiente:

######Algoritmo EM#######
fact = rep(0,100)
  for (i in 1:100){
     for (j in 1:i){
    fact[i] = fact[i] + log(j)
  }
}

retoparada = FALSE

algoritmoEM <- function (pi,data,lambds,rs,k){

while(parada == FALSE){

###### Función de densidad ########  
pmf <- function(x, lambda,r){
  lambda*exp((r-1)*log(lambda*x) - fact[r-1] - lambda*x )
}

####### E-STEP ########

#Matriz que contiene las funciones de densidad para todo m y todo k#
matrizp = matrix(0,nrow = length(data),ncol = k) 
for (i in 1:k){
  for(j in 1:length(data)){
    matrizp[j,i]= pmf(data[j],lambds[i],rs[i])
  }
}

matriznumerador = matrix(0,nrow = length(data),ncol = k) 
for (i in 1:k){
  for(j in 1:length(data)){
    matriznumerador[j,i]= (pi[i]*matrizp[j,i])
  }
}

matrizdenominador = matrix(0,nrow = length(data),ncol = 1) 
for(j in 1:length(data)){
  for (i in 1:k){
    matrizdenominador[j,1]= matrizdenominador[j,1]+matriznumerador[j,i]        
  }
}

matrizq = matrix(0,nrow = length(data),ncol = k)
for (i in 1:k){
  for(j in 1:length(data)){
    matrizq[j,i]= matriznumerador[j,i]/matrizdenominador[j,1]
  }
}

##### M-STEP ######

matriz.alfas = matrix(0, nrow = k)
matriz.lambdas = matrix(0, nrow = k)

mult <- rep(0,k)
mat.q <- rep(0,k)

for (i in 1:k)
{
  sumilla = 0
  conta = 0
for (j in 1:length(data))
{
  sumilla = sumilla + matrizq[j,i]*data[j]
  conta = conta + matrizq[j,i]
}
  mult[i] = sumilla
  mat.q[i] = conta
}

for (i in 1:k){
  matriz.alfas[i] = (1/length(data))*mat.q[i]
  matriz.lambdas[i] = (rs[i]*mat.q[i])/mult[i]
}

###### Cáculo del error ######
error <- integer
for (i in 1:k){
  error = error + (matriz.alfas[i]-pi[i])^2 + (matriz.lambdas[i]-lambds[i])^2
}
###### Condición de Parada #####
if(error <= 1e-5) {
  parada == TRUE
}

###### Actualización de Alfas y Lambdas #######
pi[i] <<- matriz.alfas[i]
lambds[i] <<- matriz.lambdas[i]
  }
}


###### DATOS ######
data1 <- (qweibull(runif(1000), shape=2.75, scale=0.25))
data1.mean <- mean(data1)
data1.mean
data1.var <- var(data1)
coefvar2 <- data1.var/(data1.mean^2)
N <- 20
  if(coefvar2<=1) {
    K1 <- 1
    K2 <- 2
    K3 <- 3
  rs1 = rep(N,K1)
  lambda1= 1/data1.mean
  pi1=1

  for(i in 1:floor(N/2)){
    for(j in floor(N/2):N-1){
      if(i+j==N){
        rs2=matrix(c(i,j),nrow=1,ncol=1)
        pi2 <- rep(1/K2,K2)

  lambda2 <- rep(0,K2)
  for(i in 1:K2)
  {
    lambda2[i] = rs2[i]/data1.mean+(1/(data1.mean*i))
  }
}
  }
 }

} else { 
   K1 <- N 
   K2 <- N-1
   K3 <- N-2
   rs1 = rep(1,N)
   rs1
   rs2 = rep(1,N-2)
   append(rs2,2,N-2)
   rs3 = rep(1,N-3)
   append(rs3,3,N-3)

   pi1 <- rep(1/K1,K1)
   pi2 <- rep(1/K2,K2)
   pi3 <- rep(1/K3,K3)

   lambda1 <- rep(0,K1)
   for(i in 1:K1)
     {
      lambda1[i] = rs1[i]/data1.mean+(1/(data1.mean*i))
     }

   lambda2 <- rep(0,K2)
   data1.mean
   for(i in 1:K2)
  {
    lambda2[i] = rs2[i]/data1.mean+(1/(data1.mean*i))
  }

  lambda3 <- rep(0,K3)
  data1.mean
  for(i in 1:K3)
  {
    lambda3[i] = rs3[i]/data1.mean+(1/(data1.mean*i))
  }
   lambda3

   variables1 <- algoritmoEM(pi1,data1,lambda1,rs1,K1)
  }

1 respuesta 1

1

El error que mencionas tiene que ver con reemplazar por ejemplo un valor en un vector con un numeric(0), en R todas las listas y vectores arrancan con índice 1, es lo que veo que te esta pasando en esta función:

pmf <- function(x, lambda, r){
      lambda*exp((r-1)*log(lambda*x) - fact[r-1] - lambda*x )
    }

Cuando usas fact[r-1], si r = 1, el resultado va a ser cero con lo cual vas a obtener un numeric(0) que se trasladara al valor de retorno.

Acá:

for (i in 1:k){
  for(j in 1:length(data)){
    matrizp[j,i]= pmf(data[j],lambds[i],rs[i])
  }
}

Justamente usando rs[i], si dicho item es 1 se producirá el error, al suceder que matrizp[j,i] = numeric(0). No sé como arreglar esta lógica por que no estoy familiarizado con el algoritmo EM.

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.