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