Skip to main content
3 de 3
Mejor explicación. Formato de código.
lois6b
  • 7.6k
  • 5
  • 33
  • 54

Error aplicando el algoritmo EM para generar una Phase Type

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)
  }