[R] Help with optim function in R, please?

[fadeh Alanazi]

Hi all, 

Many thank in advance for helping me.  I tried to fit Expectation Maximization algorithm for mixture data. I must used one of numerical method to maximize my function. 

I built my code but I do not know how to make the optim function run over a different value of the parameters.  That is, 

For E-step I need to get the value of mixture weights based on the current (initial) values of the parameter of the density. Then, multiple the weight by the logliklihood function and maximize it (M-step)
Then, I would like to take the new values of the parameter (from M-step) and plug it in the weight, to get a new value of the weight. Then, iterate till converges. 

I tried the following code, but it does not work. 

library(copula)
library(VineCopula)

## to generate the data
set.seed(123)
cp <- mixCopula(list(frankCopula(4),claytonCopula(2)))
cop <- rCopula(100,cp)
x <- pobs(cop) ## this is the data
## my function including optim function
myfunc <- function(data,copula=list(frankCopula(4,dim=2), claytonCopula(0.5,dim=2)),maxit=200){
  
  # copula[[1]]@parameters <- par[1]
  # copula[[2]]@parameters <- par[2]
  optim_1 <- function(par, data.=data, copula.=copula){
    copula[[1]]@parameters <- par[1]
    copula[[2]]@parameters <- par[2]
    tau_1 <- par[3]*dCopula(data,copula[[1]],log = F)
    tau_2 <- par[4]*dCopula(data,copula[[2]],log=F)
    Tau_1 <- tau_1 / sum(tau_1 + tau_2)
    Tau_2 <- tau_2 / sum(tau_1 + tau_2)
    up_pi1 <- sum(Tau_1)/ 100
    up_pi2 <- sum(Tau_2) / 100
    # Tau <- c(Tau_1, Tau_2)
    ll <- (-sum(Tau_1*dCopula(data,copula[[1]],log=T)))
    #ll_1 <- -sum(Tau_2*dCopula(data,copula[[2]],log=T))
    
   
    if (is.finite(ll)) {
      return(ll)
    } else {
      if (is.na(ll)) {
        message(par)
      }
      message(ll)
      return(-3^305)
    }
  }
  
  
  
  xy <- optim(par=c(2,0.5,0.2,0.2),fn=optim_1,method="L-BFGS-B",control = list(maxit= 200, trace=1),lower= c(copula[[1]]@param.lowbnd, copula[[2]]@param.lowbnd,0,0), upper = c(copula[[1]]@param.upbnd, copula[[2]]@param.upbnd,1,1))
  
  return(xy$par)
}

Any help please?

Regards, 
Fadhah 



	[[alternative HTML version deleted]]