Re: [R] Fwd: Find maxima of a function

[Martin Maechler]

> Ranjan Maitra 
> on Sun, 27 Aug 2017 08:17:15 -0500 writes:

> I have not followed the history of this thread, but I am
> quite flummoxed as to why the OP is rewriting code to
> estimate parameters from an univariate Gaussian mixture
> model when alternatives such as EMCluster (which generally
> appears to handle initialization better than MClust)
> exist. Or perhaps there is more to it in which case I
> apologize. But I thought that I would make the OP aware of
> the package (of which, in full disclosure, I am
> co-author).  Best wishes, Ranjan

in this case, I should also mention the small (but nice) package
'nor1mix' (normal 1-dimensional mixtures) for which I had added
MLE via optim() in addition to the traditional slow and
sometimes hard-to-get-to-converge-correctly EM algorithm.
After
install.packages("nor1mix")
require("nor1mix")

?norMixMLE

Martin Maechler,
ETH Zurich


> On Sun, 27 Aug 2017 16:01:59 +0300 Ismail SEZEN
>  wrote:

>> Dear Niharika,
>> 
>> As I said before, the problem is basically an
>> optimization issue. You should isolate the problematic
>> part from the rest of your study. Sometimes, more
>> information does not help to solution. All the answers
>> from us (Ulrik, David, me) are more or less are correct
>> to find a maximum point. Newton’s method is also
>> correct. But after answers, you only say, it didn’t find
>> the right maxima. At this point I’m cluesless, because
>> problem might be originated at some other point and we
>> don’t know it. For instance, In my previous answer, my
>> suggested solution was right for both your 2
>> situations. But suddenly you just said, it didin’t work
>> for a mean greater than 6 and I’m not able to reproduce
>> your new situation.
>> 
>> I’m really upset that you lost 4 weeks on a very easy
>> issue (it is very long time). But if you don’t want to
>> loose anymore, please, isolate your question from the
>> rest of the work, create minimal reproduciple example,
>> and let’s focus on the problematic part. I assure you
>> that your problem will be solved faster.
>> 
>> I could install the distr package at the end, but I’m
>> getting another error while running your code. But I
>> don’t believe the question is related to this package.
>> 
>> Let me explain about extremum points although most of you
>> know about it. Let’s assume we have a y(x) function. An
>> extremum (max/min) point is defined as where dy/dx =
>> 0. If second order derivative of y is less than 0, it’s a
>> maximum, if greater than 0, it’s a minimum point.  If
>> zero, it’s undefined. So, at the end, we need x and y
>> vectors to find maxima.
>> 
>> y <- c(1,2,3,4,3,2,3,4,5,6,7,8,9,8,7,6,5,6,7,6,5) #
>> sample y values x <- 1:length(y) # correspoinding x
>> values
>> 
>> # Now we need to convert this discrete x-y vectors to a
>> function. Because we don’t know the values between
>> points.  # That’s why, we fit a cubic spline to have the
>> values between discrete points.
>> 
>> fun <- splinefun(x = x, y = y, method = "n”)
>> 
>> # now, we are able to find what value of y at x = 1.5 by
>> the help of fun function.
>> 
>> x2 <- seq(1, max(x), 0.1) y2 <- fun(x2)
>> 
>> plot(x, y, type = "l") lines(x2, y2, col = "red”)
>> 
>> # see red and black lines are slightly different because
>> we fit a cubic spline to the discrete points.  # Now, we
>> will use optimize function to find the maxima in a
>> defined interval. Interval definition is crucial.  #
>> optimize function calculates all required steps explained
>> above to find an extremum and gives the result.
>> 
>> max.x <- optimize(fun, interval = range(x), maximum =
>> TRUE) max.x2 <- x[which(y == max(y))] # global maximum of
>> discrete x-y vectors
>> 
>> print(max.x) # x value of global maximum of y
>> print(max.x2) # x value of global maximum of y ( discrete
>> points)
>> 
>> see max.x and max.x2 are very close. max.x is calculated
>> under a cubic spline assumption and max.x2 is calculated
>> by discrete points.
>> 
>> As a result, problem is reduced to find maxima of x-y
>> vector pairs and this SHOULD work for all sort of
>> situations (UNIVERSAL FACT). If it does not work one of
>> your situations, MOST PROBABLY, you are doing something
>> wrong and we need to investigate your failed case.
>> 
>> At this point, you can use “dput” function to extract,
>> copy and paste your x-y vectors. So, we can copy the
>> statement, run in our environment without requiring distr
>> or any other package.
>> 
>> Best wishes,
>> 
>> isezen
>> 
>> 
>> > On 27 Aug 2017, at 

Re: [R] Fwd: Find maxima of a function

[Bert Gunter]

Just a quick note:

"I have tried to use optimr but it gives me the local maxima, now I am
struck with the problem of how to get the global maxima"

This is a profound misunderstanding of how numerical optimization
works. **All** numerical optimizers can provide only "local" maxima
(the exact meaning of "local" may vary). It is mathematically
impossible to numerically obtain "global" maxima and know that you
have done so.*

-- Bert

*Unless, of course, you can prove analytically that only global maxima
exist or some other such supplementary mathematical information.
Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Sat, Aug 26, 2017 at 6:39 AM, niharika singhal
 wrote:
> Hi,
>
> Thanks for your mail, and time
>
> It is not working for some arguments, when mean value is like >6.
>
>
> case
>
> mc0 <- c(0.0886,0.1744455,0.1379778,0.1209769,0.1573065,0.
> 1134463,0.2074027)
>
> rv <-UnivarMixingDistribution(Norm(486.4255, 53.24133),
>
>   Norm(664.0713, 3.674773),
>
>   Norm(669.0484, 4.101381),
>
>   Norm(677.1753, 4.869985),
>
>   Norm(683.2635, 7.288175),
>
>   Norm(727.6229, 37.64198),
>
>   Norm(819.2011, 57.06655),
>
>   mixCoeff=mc0/sum(mc0))
>
> plot(rv, to.draw.arg="d")
>
>
> I am getting 731.1345 from the code you have provide
>
>
> It is part of a code, so it was difficult to write a reproducible code
>
> I have tried to use optimr but it gives me the local maxima, now I am
> struck with the problem of how to get the global maxima
>
>
> On Sat, Aug 26, 2017 at 3:12 PM, Ismail SEZEN [via R] <
> ml+s789695n4744993...@n4.nabble.com> wrote:
>
>>
>> > On 26 Aug 2017, at 14:18, Ulrik Stervbo <[hidden email]
>> > wrote:
>> >
>> > Please keep the list in cc.
>> >
>> > Sorry, it didn't work as expected. Maybe someone else have an
>> appropriate
>> > solution.
>> >
>> > Best,
>> > Ulrik
>> >
>> > On Sa., 26. Aug. 2017, 12:57 niharika singhal <[hidden email]
>> >
>> > wrote:
>> >
>> >> Hi
>> >>
>> >> Thanks for you mail,
>> >> I really appreciate your time on my problem
>> >>
>> >> I have posted this problem on
>> >>
>> >>
>> >> https://stats.stackexchange.com/questions/299590/to-find-max
>> ima-for-gaussian-mixture-model
>> >>
>> >>
>> >> The plot I am getting using UnivarMixingDistribution from distr package
>> in
>> >> R
>> >>
>> >> code is
>> >>
>> >> mc0= c(0.1241933, 0.6329082 <06329%20082>, 0.2428986 <02428%20986>)
>> >> 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")
>> >>
>> >> I want output around 672 in first case and in 2nd case around 2.1
>> >> according to the plot.
>> >> your code will not work in both the scenario
>> >>
>> >> Regards
>> >> Niharika Singhal
>> >>
>> >>
>> >> On Sat, Aug 26, 2017 at 12:47 PM, Ulrik Stervbo <[hidden email]
>> >
>> >> wrote:
>> >>
>> >>> Hi,
>> >>>
>> >>> I once found this somewhere on stackoverflow:
>> >>>
>> >>> values <- rnorm(20, mean = c(2.15,2.0,2.9), sd = c(0.1,0.1,0.1))
>> >>>
>> >>> v_dens <- density(values)
>> >>> v_dens_y <- v_dens$y
>> >>>
>> >>> r <- rle(v_dens_y)
>> >>> # These functions ignore the extremes if they are the first or last
>> point
>> >>> maxima_index <- which(rep(x = diff(sign(diff(c(-Inf, r$values,
>> -Inf
>> >>> == -2,  times = r$lengths))
>> >>> minima_index <- which(rep(x = diff(sign(diff(c(-Inf, r$values,
>> -Inf
>> >>> == 2,  times = r$lengths))
>> >>>
>> >>> plot(v_dens_y)
>> >>>
>> >>> HTH
>> >>> Ulrik
>> >>>
>> >>>
>> >>> On Sat, 26 Aug 2017 at 11:49 niharika singhal <
>> >>> [hidden email] >
>> wrote:
>> >>>
>>  I have a Gaussian mixture model with some parameters
>> 
>>  mean=(506.8644,672.8448,829.902)
>> 
>>  sigma=(61.02859,9.149168,74.84682)
>> 
>>  c=(0.1241933, 0.6329082 <06329%20082>, 0.2428986 <02428%20986>)
>> 
>>  And the plot look something like below.[image: enter image
>> description
>>  here]
>>  
>> 
>>  Also, if I change my parameters to
>> 
>>  mean=(2.15,2.0,2.9)
>> 
>>  sigma=(0.1,0.1,0.1)
>> 
>>  c=(1/3,1/3,1/3)
>> 
>>  Then plot would change to[image: enter image description here]
>>  
>> 
>>  Is there any way to find the maxima. I have tried Newton's method but
>> it
>>  gave me the wrong output.
>> 
>>  Like in general some common solution, 

Re: [R] Fwd: Find maxima of a function

[Ranjan Maitra]

I have not followed the history of this thread, but I am quite flummoxed as to 
why the OP is rewriting code to estimate parameters from an univariate Gaussian 
mixture model when alternatives such as EMCluster (which generally appears to 
handle initialization better than MClust) exist. Or perhaps there is more to it 
in which case I apologize. But I thought that I would make the OP aware of the 
package (of which, in full disclosure, I am co-author).

Best wishes,
Ranjan


On Sun, 27 Aug 2017 16:01:59 +0300 Ismail SEZEN  wrote:

> Dear Niharika,
> 
> As I said before, the problem is basically an optimization issue. You should 
> isolate the problematic part from the rest of your study. Sometimes, more 
> information does not help to solution. All the answers from us (Ulrik, David, 
> me) are more or less are correct to find a maximum point. Newton’s method is 
> also correct. But after answers, you only say, it didn’t find the right 
> maxima. At this point I’m cluesless, because problem might be originated at 
> some other point and we don’t know it. For instance, In my previous answer, 
> my suggested solution was right for both your 2 situations. But suddenly you 
> just said, it didin’t work for a mean greater than 6 and I’m not able to 
> reproduce your new situation.
> 
>  I’m really upset that you lost 4 weeks on a very easy issue (it is very long 
> time). But if you don’t want to loose anymore, please, isolate your question 
> from the rest of the work, create minimal reproduciple example, and let’s 
> focus on the problematic part. I assure you that your problem will be solved 
> faster.
> 
> I could install the distr package at the end, but I’m getting another error 
> while running your code. But I don’t believe the question is related to this 
> package.
> 
> Let me explain about extremum points although most of you know about it. 
> Let’s assume we have a y(x) function. An extremum (max/min) point is defined 
> as where dy/dx = 0. If second order derivative of y is less than 0, it’s a 
> maximum, if greater than 0, it’s a minimum point.  If zero, it’s undefined. 
> So, at the end, we need x and y vectors to find maxima. 
> 
> y <- c(1,2,3,4,3,2,3,4,5,6,7,8,9,8,7,6,5,6,7,6,5) # sample y values
> x <- 1:length(y) # correspoinding x values
> 
> # Now we need to convert this discrete x-y vectors to a function. Because we 
> don’t know the values between points.
> # That’s why, we fit a cubic spline to have the values between discrete 
> points.
> 
> fun <- splinefun(x = x, y = y, method = "n”)
> 
> # now, we are able to find what value of y at x = 1.5 by the help of fun 
> function.
> 
> x2 <- seq(1, max(x), 0.1)
> y2 <- fun(x2)
> 
> plot(x, y, type = "l")
> lines(x2, y2, col = "red”)
> 
> # see red and black lines are slightly different because we fit a cubic 
> spline to the discrete points.
> # Now, we will use optimize function to find the maxima in a defined 
> interval. Interval definition is crucial.
> # optimize function calculates all required steps explained above to find an 
> extremum and gives the result.
> 
> max.x <- optimize(fun, interval = range(x), maximum = TRUE)
> max.x2 <- x[which(y == max(y))] # global maximum of discrete x-y vectors
> 
> print(max.x) # x value of global maximum of y
> print(max.x2) # x value of global maximum of y ( discrete points)
> 
> see max.x and max.x2 are very close. max.x is calculated under a cubic spline 
> assumption and max.x2 is calculated by discrete points.
> 
> As a result, problem is reduced to find maxima of x-y vector pairs and this 
> SHOULD work for all sort of situations (UNIVERSAL FACT). If it does not work 
> one of your situations, MOST PROBABLY, you are doing something wrong and we 
> need to investigate your failed case. 
> 
> At this point, you can use “dput” function to extract, copy and paste your 
> x-y vectors. So, we can copy the statement, run in our environment without 
> requiring distr or any other package.
> 
> Best wishes,
> 
> isezen
> 
> 
> > On 27 Aug 2017, at 12:59, niharika singhal  
> > wrote:
> > 
> > Dear David and Ismail,
> > 
> > The actual problem is I am getting the parameters from the Kmeans cluster
> > on the data set obtained from the mclust package.
> > 
> > In mclust method I am changing the value of G according to user input, so
> > the number of means, sigma and the coefficien mixture I will get is not
> > fixed, It can be 3 or 4or5 or 10. It will all depend on the user.
> > 
> > Then on the result of kmeans, I wanted to find the maxima so that I can use
> > Normalized probability formula further for my logic and in order to do that
> > I used Newton's method and that was implemented in the optimr package.
> > 
> > I was getting the wrong result so I used distr package in order to check
> > the problem and I figured out I was getting the wrong maxima. So I can to
> > the conclusion which I have posted as my query.
> > 
> > PS: I am able to use 

Re: [R] Fwd: Find maxima of a function

[Ismail SEZEN]

Dear Niharika,

As I said before, the problem is basically an optimization issue. You should 
isolate the problematic part from the rest of your study. Sometimes, more 
information does not help to solution. All the answers from us (Ulrik, David, 
me) are more or less are correct to find a maximum point. Newton’s method is 
also correct. But after answers, you only say, it didn’t find the right maxima. 
At this point I’m cluesless, because problem might be originated at some other 
point and we don’t know it. For instance, In my previous answer, my suggested 
solution was right for both your 2 situations. But suddenly you just said, it 
didin’t work for a mean greater than 6 and I’m not able to reproduce your new 
situation.

 I’m really upset that you lost 4 weeks on a very easy issue (it is very long 
time). But if you don’t want to loose anymore, please, isolate your question 
from the rest of the work, create minimal reproduciple example, and let’s focus 
on the problematic part. I assure you that your problem will be solved faster.

I could install the distr package at the end, but I’m getting another error 
while running your code. But I don’t believe the question is related to this 
package.

Let me explain about extremum points although most of you know about it. Let’s 
assume we have a y(x) function. An extremum (max/min) point is defined as where 
dy/dx = 0. If second order derivative of y is less than 0, it’s a maximum, if 
greater than 0, it’s a minimum point.  If zero, it’s undefined. So, at the end, 
we need x and y vectors to find maxima. 

y <- c(1,2,3,4,3,2,3,4,5,6,7,8,9,8,7,6,5,6,7,6,5) # sample y values
x <- 1:length(y) # correspoinding x values

# Now we need to convert this discrete x-y vectors to a function. Because we 
don’t know the values between points.
# That’s why, we fit a cubic spline to have the values between discrete points.

fun <- splinefun(x = x, y = y, method = "n”)

# now, we are able to find what value of y at x = 1.5 by the help of fun 
function.

x2 <- seq(1, max(x), 0.1)
y2 <- fun(x2)

plot(x, y, type = "l")
lines(x2, y2, col = "red”)

# see red and black lines are slightly different because we fit a cubic spline 
to the discrete points.
# Now, we will use optimize function to find the maxima in a defined interval. 
Interval definition is crucial.
# optimize function calculates all required steps explained above to find an 
extremum and gives the result.

max.x <- optimize(fun, interval = range(x), maximum = TRUE)
max.x2 <- x[which(y == max(y))] # global maximum of discrete x-y vectors

print(max.x) # x value of global maximum of y
print(max.x2) # x value of global maximum of y ( discrete points)

see max.x and max.x2 are very close. max.x is calculated under a cubic spline 
assumption and max.x2 is calculated by discrete points.

As a result, problem is reduced to find maxima of x-y vector pairs and this 
SHOULD work for all sort of situations (UNIVERSAL FACT). If it does not work 
one of your situations, MOST PROBABLY, you are doing something wrong and we 
need to investigate your failed case. 

At this point, you can use “dput” function to extract, copy and paste your x-y 
vectors. So, we can copy the statement, run in our environment without 
requiring distr or any other package.

Best wishes,

isezen


> On 27 Aug 2017, at 12:59, niharika singhal  
> wrote:
> 
> Dear David and Ismail,
> 
> The actual problem is I am getting the parameters from the Kmeans cluster
> on the data set obtained from the mclust package.
> 
> In mclust method I am changing the value of G according to user input, so
> the number of means, sigma and the coefficien mixture I will get is not
> fixed, It can be 3 or 4or5 or 10. It will all depend on the user.
> 
> Then on the result of kmeans, I wanted to find the maxima so that I can use
> Normalized probability formula further for my logic and in order to do that
> I used Newton's method and that was implemented in the optimr package.
> 
> I was getting the wrong result so I used distr package in order to check
> the problem and I figured out I was getting the wrong maxima. So I can to
> the conclusion which I have posted as my query.
> 
> PS: I am able to use distr package without any problem.
> 
> And I want to use dnorm because I have a Gaussian mixture model, so I did
> not want to use rnorm.
> 
> I hope you get what my problem is now.
> 
> I will try both of your solution, which uses the same function and sees if
> it helps me.
> 
> I am struck at this point from almost 4 weeks and it is delaying my work a
> lot.
> 
> I really appreciate all your help
> 
> Thanks & Regards
> Niharika SInghal

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