[R] Fwd: Find maxima of a function
Martin Maechler
maechler at stat.math.ethz.ch
Tue Aug 29 14:59:10 CEST 2017
>>>>> Ranjan Maitra <maitra at email.com>
>>>>> 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
> <sezenismail at gmail.com> 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
>> <niharikasinghal1990 at gmail.com> 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 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.
> --
> Important Notice: This mailbox is ignored: e-mails are set
> to be deleted on receipt. Please respond to the mailing
> list if appropriate. For those needing to send personal or
> professional e-mail, please use appropriate addresses.
> ______________________________________________
> 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.
More information about the R-help
mailing list