Luigi Marongiu
2014-Apr-21 18:24 UTC
[R] find maximum values on the density function of a bimodal distribution
*dear all,I have a bimodal distribution and i would like to identify the two most frequent values. Using unimodal values i found from the R archive that is possible to identify the most frequent value, which corresponds to the peak of the density function:from Bert Gunter, Fri Mar 13 04rv <- rbinom(10000,1,0.1) + rnorm(10000)d.rv = density(rv)d.x = d.rv$xd.y d.rv$yd.rv.max = d.rv$x[which.max(d.rv$y)]plot(d.rv)abline(v=d.rv.max)The following is instead a bimodal distribution, so how to identify the two peak values?rv <-c(-0.161948986691092, 6.25551346546337, 14.4449902518518, 5.71035851092391, 14.2659690175213, 13.4326694168529, 2.25851352934142, 19.0283322756163, 18.7578638985237, 3.6079681449722, 12.8469831785319, 16.8172446560292, 0.586180794938773, -2.78468096182699, 6.41435801855101, 5.40404032358274, 16.422996950473, 0.0146677573458032, 6.42413783291763, 16.9133720373153, 21.2106745914211, 4.04126872437582, 10.8506381906698, 8.74969277743266, 14.2820697977719, 15.9914414284944, 8.35395871093583, 22.322063211793, 14.3922587603495, -0.889640152783791, 15.5590006991235, 13.8636215566311, 8.04175502056292, -1.85932938527655, 3.0791620072771, 20.5240745935955, 3.68821147789088, 7.38116287748327, 8.13585591202069, 5.94733543506892, 7.77891267272758, 14.4774347329043, 0.309628817532129, -9.0260821589798, 12.4102239527495, -4.64595423395751, 17.3141235128072, 20.1952807795295, -8.18424754421263, -1.58917260547841)d.rv density(rv)d.x = d.rv$xd.y = d.rv$yd.rv.max d.rv$x[which.max(d.rv$y)]plot(d.rv)Thank you.Best wishes,Luigi* [[alternative HTML version deleted]]
David Carlson
2014-Apr-21 20:38 UTC
[R] find maximum values on the density function of a bimodal distribution
This will work, as long as there are exactly 2 distinct modes:> runs <- rle(sign(diff(d.rv$y))) > length(runs$lengths) # There should be 4 runs if there areexactly 2 modes [1] 4> mode1 <- runs$lengths[1]+1 > mode2 <- length(d.rv$x)- runs$lengths[4] > d.rv$y[c(mode1, mode2)] # The 2 modes:[1] 0.04012983 0.04049404 -------------------------------------------------- David L Carlson Listowner, Webmaster http://people.tamu.edu/~dcarlson/nar/ ----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Luigi Marongiu Sent: Monday, April 21, 2014 1:24 PM To: r-help at r-project.org Subject: [R] find maximum values on the density function of a bimodal distribution *dear all,I have a bimodal distribution and i would like to identify the two most frequent values. Using unimodal values i found from the R archive that is possible to identify the most frequent value, which corresponds to the peak of the density function:from Bert Gunter, Fri Mar 13 04rv <- rbinom(10000,1,0.1) + rnorm(10000)d.rv = density(rv)d.x d.rv$xd.y d.rv$yd.rv.max d.rv$x[which.max(d.rv$y)]plot(d.rv)abline(v=d.rv.max)The following is instead a bimodal distribution, so how to identify the two peak values?rv <-c(-0.161948986691092, 6.25551346546337, 14.4449902518518, 5.71035851092391, 14.2659690175213, 13.4326694168529, 2.25851352934142, 19.0283322756163, 18.7578638985237, 3.6079681449722, 12.8469831785319, 16.8172446560292, 0.586180794938773, -2.78468096182699, 6.41435801855101, 5.40404032358274, 16.422996950473, 0.0146677573458032, 6.42413783291763, 16.9133720373153, 21.2106745914211, 4.04126872437582, 10.8506381906698, 8.74969277743266, 14.2820697977719, 15.9914414284944, 8.35395871093583, 22.322063211793, 14.3922587603495, -0.889640152783791, 15.5590006991235, 13.8636215566311, 8.04175502056292, -1.85932938527655, 3.0791620072771, 20.5240745935955, 3.68821147789088, 7.38116287748327, 8.13585591202069, 5.94733543506892, 7.77891267272758, 14.4774347329043, 0.309628817532129, -9.0260821589798, 12.4102239527495, -4.64595423395751, 17.3141235128072, 20.1952807795295, -8.18424754421263, -1.58917260547841)d.rv density(rv)d.x = d.rv$xd.y = d.rv$yd.rv.max d.rv$x[which.max(d.rv$y)]plot(d.rv)Thank you.Best wishes,Luigi* [[alternative HTML version deleted]] ______________________________________________ R-help at r-project.org mailing list 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.