George W. Gilchrist
2004-Sep-16 14:59 UTC
[R] Estimating parameters for a bimodal distribution
For several years, I have been using Splus to analyze an ongoing series of
datasets that have a bimodal distribution. I have used the following
functions, in particular the ms() function, to estimate the parameters: two
means, two standard deviations, and one proportion. Here is the code I've
been using in S:
btmp.bi <- function(vec, p, m1, m2, sd1, sd2)
{
(exp(p)/(1+exp(p)))*dnorm(vec,mean=m1,sd=abs(sd1))+
(1-(exp(p)/(1+exp(p))))*dnorm(vec,mean=m2,sd=abs(sd2))
}
btmp11 <- ms( ~ - sum(log((btmp.bi(btmp1$Temp, p, m1, m2,
s1, s2)))), start = list(p = 0.4, m1 = 38, m2 = 40, s1
= 1, s2 = 1), control = list(maxiter = 200))
I have looked in the archives and tried various alternatives, especially
optim(), but so far have had nothing but frustration. I've been running this
in a semi-automated program on several hundred datasets a few times each
month for several years now. I would like to figure out how to move this to
R. Thank you for any help you might offer.
=================================================================George W.
Gilchrist Email #1: gwgilc at wm.edu
Department of Biology, Box 8795 Email #2: kitesci at cox.net
College of William & Mary Phone: (757) 221-7751
Williamsburg, VA 23187-8795 Fax: (757) 221-6483
http://gwgilc.people.wm.edu/
Dimitris Rizopoulos
2004-Sep-16 15:42 UTC
[R] Estimating parameters for a bimodal distribution
Hi George,
I tried the following and it worked for me,
btmp.bi <- function(par., vec){
p <- par.[1]
mu1 <- par.[2]
mu2 <- par.[3]
sigma1 <- par.[4]
sigma2 <- par.[5]
-sum(log( plogis(p)*dnorm(vec, mu1, abs(sigma1)) +
(1-plogis(p))*dnorm(vec, mu2, abs(sigma2)) ))
}
vec <- ifelse(runif(500)<.4, rnorm(1, 38, 1), rnorm(1, 40, 1))
optim(par=c(p=0.4, mu1=38, mu2=40, sigma1=1, sigma2=1), btmp.bi,
method="CG", vec=vec)
I hope it helps.
Best,
Dimitris
----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/396887
Fax: +32/16/337015
Web: http://www.med.kuleuven.ac.be/biostat/
http://www.student.kuleuven.ac.be/~m0390867/dimitris.htm
----- Original Message -----
From: "George W. Gilchrist" <gwgilc at wm.edu>
To: <r-help at stat.math.ethz.ch>
Sent: Thursday, September 16, 2004 4:59 PM
Subject: [R] Estimating parameters for a bimodal distribution
> For several years, I have been using Splus to analyze an ongoing
> series of
> datasets that have a bimodal distribution. I have used the following
> functions, in particular the ms() function, to estimate the
> parameters: two
> means, two standard deviations, and one proportion. Here is the code
> I've
> been using in S:
>
> btmp.bi <- function(vec, p, m1, m2, sd1, sd2)
> {
> (exp(p)/(1+exp(p)))*dnorm(vec,mean=m1,sd=abs(sd1))+
> (1-(exp(p)/(1+exp(p))))*dnorm(vec,mean=m2,sd=abs(sd2))
> }
> btmp11 <- ms( ~ - sum(log((btmp.bi(btmp1$Temp, p, m1, m2,
> s1, s2)))), start = list(p = 0.4, m1 = 38, m2 = 40, s1
> = 1, s2 = 1), control = list(maxiter = 200))
>
> I have looked in the archives and tried various alternatives,
> especially
> optim(), but so far have had nothing but frustration. I've been
> running this
> in a semi-automated program on several hundred datasets a few times
> each
> month for several years now. I would like to figure out how to move
> this to
> R. Thank you for any help you might offer.
>
> =================================================================>
George W. Gilchrist Email #1: gwgilc at wm.edu
> Department of Biology, Box 8795 Email #2: kitesci at cox.net
> College of William & Mary Phone: (757) 221-7751
> Williamsburg, VA 23187-8795 Fax: (757) 221-6483
> http://gwgilc.people.wm.edu/
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
George,
Venables & Ripley: Modern Applied Statistics with S, Springer 2002, Chapter
16 (An example: fitting a mixture model) may be helpful.
Hannu Kahra
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch]On Behalf Of George W.
Gilchrist
Sent: Thursday, September 16, 2004 5:00 PM
To: r-help at stat.math.ethz.ch
Subject: [R] Estimating parameters for a bimodal distribution
For several years, I have been using Splus to analyze an ongoing series of
datasets that have a bimodal distribution. I have used the following
functions, in particular the ms() function, to estimate the parameters: two
means, two standard deviations, and one proportion. Here is the code I've
been using in S:
btmp.bi <- function(vec, p, m1, m2, sd1, sd2)
{
(exp(p)/(1+exp(p)))*dnorm(vec,mean=m1,sd=abs(sd1))+
(1-(exp(p)/(1+exp(p))))*dnorm(vec,mean=m2,sd=abs(sd2))
}
btmp11 <- ms( ~ - sum(log((btmp.bi(btmp1$Temp, p, m1, m2,
s1, s2)))), start = list(p = 0.4, m1 = 38, m2 = 40, s1
= 1, s2 = 1), control = list(maxiter = 200))
I have looked in the archives and tried various alternatives, especially
optim(), but so far have had nothing but frustration. I've been running this
in a semi-automated program on several hundred datasets a few times each
month for several years now. I would like to figure out how to move this to
R. Thank you for any help you might offer.
=================================================================George W.
Gilchrist Email #1: gwgilc at wm.edu
Department of Biology, Box 8795 Email #2: kitesci at cox.net
College of William & Mary Phone: (757) 221-7751
Williamsburg, VA 23187-8795 Fax: (757) 221-6483
http://gwgilc.people.wm.edu/
______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Prof Brian Ripley
2004-Sep-16 15:57 UTC
[R] Estimating parameters for a bimodal distribution
MASS4, chapter 16 has examples of a mixture of two normals by various methods, very like this one, in both S-PLUS and R. On Thu, 16 Sep 2004, George W. Gilchrist wrote:> For several years, I have been using Splus to analyze an ongoing series of > datasets that have a bimodal distribution. I have used the following > functions, in particular the ms() function, to estimate the parameters: two > means, two standard deviations, and one proportion. Here is the code I've > been using in S: > > btmp.bi <- function(vec, p, m1, m2, sd1, sd2) > { > (exp(p)/(1+exp(p)))*dnorm(vec,mean=m1,sd=abs(sd1))+ > (1-(exp(p)/(1+exp(p))))*dnorm(vec,mean=m2,sd=abs(sd2)) > } > btmp11 <- ms( ~ - sum(log((btmp.bi(btmp1$Temp, p, m1, m2, > s1, s2)))), start = list(p = 0.4, m1 = 38, m2 = 40, s1 > = 1, s2 = 1), control = list(maxiter = 200)) > > I have looked in the archives and tried various alternatives, especially > optim(), but so far have had nothing but frustration. I've been running this > in a semi-automated program on several hundred datasets a few times each > month for several years now. I would like to figure out how to move this to > R. Thank you for any help you might offer. > > =================================================================> George W. Gilchrist Email #1: gwgilc at wm.edu > Department of Biology, Box 8795 Email #2: kitesci at cox.net > College of William & Mary Phone: (757) 221-7751 > Williamsburg, VA 23187-8795 Fax: (757) 221-6483 > http://gwgilc.people.wm.edu/ > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html > >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595