[R] simple mixture
Bill Simpson
wsi at gcal.ac.uk
Thu Nov 9 17:41:41 CET 2000
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list