Error applying the EM algorithm to generate a Phase Type

From some data I am generating (weibull) I am trying to make a Phase Type fitting using the MS algorithm in R. However, when generating the density function (matrixp) I get the following error:

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

The code is as follows:

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