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

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

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

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

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

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

"matrizp[j,i]= pmf(data[j],lambds[i],rs[i]) replacement has length zero"El código es el siguiente:

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

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

El código es el siguiente:

Origen Enlace

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"

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