R help - Aug 2017 - Problem in optimization of Gaussian Mixture model

Hello,


I am facing a problem with optimization in R from 2-3 weeks.



I have some Gaussian mixtures parameters and I want to find the maximum in
that



*Parameters are in the form *

mean1    mean2    mean3   sigma1   sigma2   sigma3   c1        c2        c3

506.8644 672.8448 829.902 61.02859 9.149168 74.84682 0.1241933
0.6329082 0.2428986





I have used optima and optimx to find the maxima, but it gives me value
near by the highest mean as an output, for example 830 in the above
parameters. The code for my optim function is

val1=*optim*(par=c(X=1000, *seg=1*),fn=xnorm,  opt=c(NA,seg=1),
method="BFGS",lower=-Inf,upper=+Inf, control=list(fnscale=-1))

*I am running the optim function in a loop and for different initial value
and taking the val1$par[1] as the best value.*

*I am taking this parameter since I want to run it on n number of segments
latter*

and xnorm is simply calculating the dnorm

*xnorm*=function(param, opt = rep(NA, 2)){

  if (any(!sapply(opt, is.na))) {

    i = !sapply(opt, is.na)

    # Fix non-NA values

    param[i] <- opt[i]

  }

  xval= param[1]

    seg <- param[2]

  sum_prob=0

  val=0

l=3

meanval=c(506.8644, 672.8448, 829.902)

sigmaval=c(61.02859, 9.149168, 74.84682)

coeffval(0.1241933, 0.6329082, 0.2428986)

  for(n in 1 :l)

  {

    mu=meanval[seg,n]

    sg=sigmaval[seg,n]

    cval=coeffval[seg,n]

    val=cval*(dnorm(xval,mu,sg))

    #print(paste0("The dnorm value for x is.: ", val))

    sum_prob=sum_prob+val

  }

  sum_prob

}


The output is not correct. Since I check my data using*
UnivarMixingDistribution* from distr package and according to this the max
should lie somewhere between 600-800

Code I used to check

mc0=c( 0.6329082,0.6329082,0.2428986)

rv
<-UnivarMixingDistribution(Norm(506.8644,61.02859),Norm(672.8448,9.149168),Norm(
829.902,74.84682), mixCoeff=mc0/sum(mc0))

plot(rv, to.draw.arg="d")



Can someone please help how I can solve this problem?


Thanks & Regards

Niharika Singhal

	[[alternative HTML version deleted]]

You are doing yourself no favors by posting HTML format email... the code 
shown below has extra characters in it that you probably did not put 
there. This is a plain text mailing list, so you need to send plain text 
email. Read the Posting Guide mentioned in the footer of this and every 
R-help posting.

On top of that, your code is not presented as a stand-alone reproducible 
example. Read [1], [2], and in particular [3]. Failing to make it 
reproducible means most readers will just ignore it and move on.

I cannot figure out why you have "seg" as an optimized parameter... in
particular, you are trying to index one-dimensional vectors with two 
indexes, one of which is the floating point variable "seg".

Anyway, despite the above, here is my take on your problem.

####
myfunc <- function( x ) {
    meanval <- c( 506.8644, 672.8448, 829.902 )
    sigmaval <- c( 61.02859, 9.149168, 74.84682 )
    coeffval <- c( 0.1241933, 0.6329082, 0.2428986 )
    sapply( x
          , function( xi )
              sum( coeffval * dnorm( xi
                                   , meanval
                                   , sigmaval
                                   )
                 )
          )
}

plot( 400:1000, myfunc( 400:1000 ) )

#' ![](http://i.imgur.com/65fhqrL.png)


# fooled by local maximum
val1 <- optim( par = c( x = 800 )
              , fn = myfunc
              , method = "BFGS"
              , control = list( fnscale= -1
                              , parscale = 1/0.025
                              )
              )
val1
#> $par
#>        x
#> 829.5249
#>
#> $value
#> [1] 0.001294662
#>
#> $counts
#> function gradient
#>        5        4
#>
#> $convergence
#> [1] 0
#>
#> $message
#> NULL

val2 <- optim( par = c( x = 1000 )
              , fn = myfunc
              , method = "SANN"
              , control = list( fnscale= -1
                              , parscale = 1/0.025
                              )
              )
val2
#> $par
#>        x
#> 672.8166
#>
#> $value
#> [1] 0.02776057
#>
#> $counts
#> function gradient
#>    10000       NA
#>
#> $convergence
#> [1] 0
#>
#> $message
#> NULL
####

[1]
http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example

[2] http://adv-r.had.co.nz/Reproducibility.html

[3] https://cran.r-project.org/web/packages/reprex/index.html (read the 
vignette)

On Thu, 24 Aug 2017, niharika singhal wrote:
> Hello,
>
>
> I am facing a problem with optimization in R from 2-3 weeks.
>
>
>
> I have some Gaussian mixtures parameters and I want to find the maximum in
> that
>
>
>
> *Parameters are in the form *
>
> mean1    mean2    mean3   sigma1   sigma2   sigma3   c1        c2        c3
>
> 506.8644 672.8448 829.902 61.02859 9.149168 74.84682 0.1241933
> 0.6329082 0.2428986
>
>
>
>
>
> I have used optima and optimx to find the maxima, but it gives me value
> near by the highest mean as an output, for example 830 in the above
> parameters. The code for my optim function is
>
> val1=*optim*(par=c(X=1000, *seg=1*),fn=xnorm,  opt=c(NA,seg=1),
> method="BFGS",lower=-Inf,upper=+Inf, control=list(fnscale=-1))
>
> *I am running the optim function in a loop and for different initial value
> and taking the val1$par[1] as the best value.*
>
> *I am taking this parameter since I want to run it on n number of segments
> latter*
>
> and xnorm is simply calculating the dnorm
>
> *xnorm*=function(param, opt = rep(NA, 2)){
>
>  if (any(!sapply(opt, is.na))) {
>
>    i = !sapply(opt, is.na)
>
>    # Fix non-NA values
>
>    param[i] <- opt[i]
>
>  }
>
>  xval= param[1]
>
>    seg <- param[2]
>
>  sum_prob=0
>
>  val=0
>
> l=3
>
> meanval=c(506.8644, 672.8448, 829.902)
>
> sigmaval=c(61.02859, 9.149168, 74.84682)
>
> coeffval(0.1241933, 0.6329082, 0.2428986)
>
>  for(n in 1 :l)
>
>  {
>
>    mu=meanval[seg,n]
>
>    sg=sigmaval[seg,n]
>
>    cval=coeffval[seg,n]
>
>    val=cval*(dnorm(xval,mu,sg))
>
>    #print(paste0("The dnorm value for x is.: ", val))
>
>    sum_prob=sum_prob+val
>
>  }
>
>  sum_prob
>
> }
>
>
> The output is not correct. Since I check my data using*
> UnivarMixingDistribution* from distr package and according to this the max
> should lie somewhere between 600-800
>
> Code I used to check
>
> mc0=c( 0.6329082,0.6329082,0.2428986)
>
> rv
>
<-UnivarMixingDistribution(Norm(506.8644,61.02859),Norm(672.8448,9.149168),Norm(
> 829.902,74.84682), mixCoeff=mc0/sum(mc0))
>
> plot(rv, to.draw.arg="d")
>
>
>
> Can someone please help how I can solve this problem?
>
>
> Thanks & Regards
>
> Niharika Singhal
>
> 	[[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.
>
---------------------------------------------------------------------------
Jeff Newmiller                        The     .....       .....  Go Live...
DCN:<jdnewmil at dcn.davis.ca.us>        Basics: ##.#.       ##.#.  Live
Go...
                                       Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
/Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k