[R] Estimating parameters for a bimodal distribution

[Kahra Hannu]

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