[R] Estimating parameters for a bimodal distribution

[George W. Gilchrist]

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/