[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 <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