I haven't tried your code, but I've seen far too many attempts like this where there is no simple call to the function with starting parameters to see if a sensible answer is returned. If you get NaN or Inf or ... that isn't sensible, then you know it is not a good idea to try further.
Beyond this, you don't have gradients. I think I would try nmkb (an updated Nelder-Mead style code that allows for bounds) from dfoptim package. Note that you cannot start on the bounds with this code. And it won't be terribly efficient, but it will -- if you ask for intermediate results -- give you some idea about your function. If everything then working reasonably, you can try L-BFGS-B, possibly with "better" starting parameters. However, I'm fairly certain you will find something to fix before then. It really does pay to be very pedestrian in using optimization software. I've been at this game (and note I wrote 3 of the 5 routines in optim(), which now all have replacements!) and I am constantly reminded in my own work to ALWAYS run some simple tests with the objective function and, if I have them, the gradients. JN On 2017-08-06 11:11 AM, fadeh Alanazi wrote:
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]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
