I presumed there was some reason why the 'optimise' function did not suit you.
optimize(function (x) a*dnorm(x, x1, s1) + (1-a)*dnorm(x, x2, s2), range(c(x1, x2))) should do the trick. On Thu, 14 Oct 2021 at 23:42, Luigi Marongiu <[email protected]> wrote: > > Thank you, I hoped there might be an automatic method more than > function analysis... > > On Thu, Oct 14, 2021 at 9:06 AM Richard O'Keefe <[email protected]> wrote: > > > > Do you really want the minimum? > > It sounds as though your model is a*N(x1,s1) + (1-a)*N(x2,s2) where > > you use mixtools to estimate > > the parameters. Finding the derivative of that is fairly > > straightforward calculus, and solving for the > > derivative being zero gives you extrema (you want the one between x1 and > > x2). > > > > In practice I might start with a quick and dirty hack like > > x <- seq(from min(x1,x2), to=max(x1,x2), length=399) > > p <- a*dnorm(x, x1, s1) + (1-a)*dnorm(x, x2, s2) > > m <- min(p) > > x[x == m] > > > > On Wed, 13 Oct 2021 at 22:12, Luigi Marongiu <[email protected]> > > wrote: > > > > > > Hello, > > > I have two peaks distributions (and some noise values in between), and > > > I used mixtools to characterize the two overlapping data. How can I > > > find the minimum between the peak distributions? > > > The sigma value of the second population is close to that value, but I > > > am not sure if this is really the cut-off point between the two > > > distributions (by eye, I would say the value is 0.125 against 0.1 of > > > the sigma value). Is there a better approach? > > > Thanks > > > > > > ``` > > > set.seed(10) > > > negative_mr <- runif(50, 0, 0.099) > > > negative_fcn <- runif(50, 1, 40) > > > positive_mr <- c(runif(30, 0.2, 0.5), runif(20, 0.4, 0.5)) > > > positive_fcn <- c(runif(30, 25, 40), runif(20, 10, 25)) > > > uncertain_mr <- runif(10, 0.099, 0.2) > > > uncertain_fcn <- runif(10, 2, 40) > > > df <- data.frame(Y=c(negative_mr, uncertain_mr, positive_mr), > > > X=c(negative_fcn, uncertain_fcn, positive_fcn), > > > GROUP=c(rep("negative", length(negative_mr)), > > > rep("uncertain", length(uncertain_mr)), > > > rep("positive", length(positive_mr)))) > > > library(mixtools) > > > model = normalmixEM((x = df$Y)) > > > summary(model) > > > # plot > > > plot(model, which=2) > > > Cut_off <- model$sigma[2] > > > Cut_offE <- 0.125 > > > abline(v=Cut_off, col="blue", lwd=2) > > > abline(v=Cut_offE, col="blue", lwd=2, lty=2) > > > plot(df$Y[df$GROUP=="negative"]~df$X[df$GROUP=="negative"], > > > pch=1, ylim=c(0,0.5), xlim=c(0,41), cex=1.5, > > > xlab=expression(bold("X")), > > > ylab=expression(bold("Y"))) > > > points(df$Y[df$GROUP=="positive"]~df$X[df$GROUP=="positive"], > > > pch=16, cex=1.5) > > > points(df$Y[df$GROUP=="uncertain"]~df$X[df$GROUP=="uncertain"], > > > pch=16, cex=1.5, col="grey") > > > abline(h=Cut_off, col="blue", lwd=2) > > > abline(h=Cut_offE, col="blue", lwd=2, lty=2) > > > ``` > > > > > > > > > -- > > > Best regards, > > > Luigi > > > > > > ______________________________________________ > > > [email protected] 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. > > > > -- > Best regards, > Luigi ______________________________________________ [email protected] 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.
