[R] simple mixture

[Bill Simpson]

On Thu, 9 Nov 2000, Emmanuel Paradis wrote:
> I am trying to do some simple mixture analyses. For instance, I have a
> sample of n observations and I suspect they come from two different
> exponential distributions with parameters rate1 and rate2, respectively.
> So, I want to estimate rate1, rate2, and the proportions of both kinds of
> individuals in the sample. I had a look at the packages mda and mclust, but
> they do not seem to do this kind of analysis.
This is an MLE problem; use nlm().

I recently had to do this with the Gumbel distribution. Experts, please
let me know if I've screwed up.

Bill

==================================================
pgumbel<-function(x,a,b) exp(-exp(-(x-a)/b))

dgumbel<-function(x,a,b) (1/b)*exp(-(x-a)/b)*exp(-exp(-(x-a)/b))

rgumbel<-function(n,a,b) a-b*log(-log(runif(n=n)))

dmixgumbel<-function(x,p,a1,b1,a2,b2)
{
p*(1/b1)*exp(-(x-a1)/b1)*exp(-exp(-(x-a1)/b1)) +
(1-p)*(1/b2)*exp(-(x-a2)/b2)*exp(-exp(-(x-a2)/b2))
}

#not sure if this is right!!!!!!
pmixgumbel<-function(x,p,a1,b1,a2,b2)
{
p*exp(-exp(-(x-a1)/b1)) + (1-p)*exp(-exp(-(x-a2)/b2))
}

fitmixture.sim<-function(p=.9,a1=200,b1=50,a2=600,b2=20)
{
n1<-round(p*150)
n2<-150-n1
x<-c(rgumbel(n=n1,a=a1,b=b1),rgumbel(n=n2,a=a2,b=b2))
plot(sort(x),ppoints(x))

fn<-function(p)
        {
        sum(-log(dmixgumbel(x,p=p[1],a1=p[2],b1=p[3],a2=p[4],b2=p[5])) )
        }
out<-nlm(fn,p=c(.9,200,50,600,20), hessian=TRUE)

xfit<-seq(min(x),max(x),length=40)
yfit<-pmixgumbel(xfit,p=out$estimate[1],a1=out$estimate[2],
b1=out$estimate[3],a2=out$estimate[4],b2=out$estimate[5])
lines(xfit,yfit)

list(out=out)
}


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._