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 at 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.
>