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: '[EMAIL PROTECTED]'; '[EMAIL PROTECTED]'
> 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
