Re: [R] finite mixture model (2-component gaussian): plotting component gaussian components?

From: <Karen.Green_at_sanofi-aventis.com>
Date: Tue 25 Oct 2005 - 01:37:53 EST


Dear R Community,

I was finally able to figure out how to plot the components of the finite mixture model.

Because it was not obvious, I am including the script here -- in case anyone in the future also needs to do this.

Karen

---
Karen M. Green, Ph.D.
Research Investigator
Drug Design Group
Sanofi Aventis Pharmaceuticals
Tucson, AZ, USA


#########################################################
# name: testMclust1Dplot
#
# date: 23 October 2005
#
# status: validated
#
# dependencies: package mclust
#
# purpose: to test the mclust1Dplot function &
# to test the plot of components
#########################################################
# (0) preliminaries
##############################################
library(mclust)
##############################################
# (1) get some data
##############################################
data(faithful) attach(faithful) eruptions<-sort(eruptions) eruptions<-data.frame(eruptions) eruptions<-as.numeric(eruptions$eruptions) detach(faithful)
##############################################
# (2) get model
##############################################
MyMixtureModel<-summary(EMclust(eruptions),eruptions)
##############################################
# (3) plot mixture model
##############################################
attach(MyMixtureModel) mclust1Dplot(data=eruptions,z=z,mu=mu,sigmasq=sigmasq,pro=pro,ask=FALSE,type=c("density")) do.call("mclust1Dplot",c(list(data=eruptions,ask=FALSE,type=c("density")),MyMixtureModel)) detach(MyMixtureModel)
##############################################
# (4) plot components
##############################################
temp <- cdens(modelName = MyMixtureModel$modelName,data=eruptions,mu=MyMixtureModel$mu, sigmasq=MyMixtureModel$sigmasq,decomp=MyMixtureModel$decomp) MyMixtureComponents<-cbind(eruptions,temp) colnames(MyMixtureComponents)<-c("eruptions","component1","component2","component3","component4") MyMixtureComponents<-data.frame(MyMixtureComponents) MyMixtureComponents$eruptions<-as.numeric(MyMixtureComponents$eruptions) attach(MyMixtureComponents) lines(x=eruptions,y=component1*MyMixtureModel$pro[1],type="l",col="red") lines(x=eruptions,y=component2*MyMixtureModel$pro[2],type="l",col="red") lines(x=eruptions,y=component3*MyMixtureModel$pro[3],type="l",col="red") lines(x=eruptions,y=component4*MyMixtureModel$pro[4],type="l",col="red") detach(MyMixtureComponents)
#########################################################
> -----Original Message----- > From: Green, Karen M. PH/US > Sent: Thursday, October 20, 2005 6:42 PM > To: 'R-help@lists.R-project.org'; 'r-help-request@stat.math.ethz.ch' > Cc: Green, Karen M. PH/US > Subject: finite mixture model (2-component gaussian): plotting component gaussian components? > Importance: High > > Dear Knowledgeable R Community Members, > > Please excuse my ignorance, I apologize in advance if this is an easy question, but I am a bit stumped and could use a little guidance. > > I have a finite mixture modeling problem -- for example, a 2-component gaussian mixture -- where the components have a large overlap, and > I am trying to use the "mclust" package to solve this problem. > > I need to decompose that mixture into its 2 components which will need to be plotted. > > What I don't know how to do is: > > (1) restrict the number of components to 2 in the "EMclust" function > > (2) obtain and plot a component gaussian density > > Regarding (1), I think this should be the 'G' value but the documentation is somewhat cryptic as regards the format. > Regarding (2), I think I need to use the "cdens" function, but again the documentation is somewhat cryptic. > > Here is a little test script to illustrate. (Note: my real dataset will not have peaks this well separated, but I needed to find a small example.)> > > ################## > # (1) get some data > ################## > data(faithful) > > ################## > # (2) get model > ################## > library(mclust) > MyMixtureModel<-summary(EMclust(faithful$eruptions),faithful$eruptions) > > ################## > # (3) plot mixture model > ################## > attach(MyMixtureModel) > mclust1Dplot(data=faithful$eruptions,z=z,mu=mu,sigmasq=sigmasq,pro=pro,ask=FALSE,type=c("density")) > do.call("mclust1Dplot",c(list(data=faithful$eruptions,ask=FALSE,type=c("density")),MyMixtureModel)) > > ################## > # (4) plot components > ################## > ??? > > Any information you might be able to shed on this would be very much appreciated. > > With appreciation for your help, > > Karen > --- > Karen M. Green, Ph.D. > Research Investigator > Drug Design Group > Sanofi Aventis Pharmaceuticals > Tucson, AZ > USA > [[alternative HTML version deleted]] ______________________________________________ R-help@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
Received on Tue Oct 25 01:45:18 2005

This archive was generated by hypermail 2.1.8 : Tue 25 Oct 2005 - 03:14:40 EST