[R] Estimating parameters for a bimodal distribution

[Dimitris Rizopoulos]

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
>