2009/6/3 Romain François <francoisrom...@free.fr>: > Hi, > > Would you like this to go to the graph gallery ? > http://addictedtor.free.fr/graphiques/ >
I would be honored to have a graph in your fine display. Do you want me to submit it formally by emailing you as described here: http://wiki.r-project.org/rwiki/doku.php?id=graph_gallery:graph_gallery I'm still campaigning to add this to ?plotmath because it is a tricky problem to combine several variable values and plotmath into a single expression. Consider the graph gallery 2d normal model http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=42 When I saw the subtitle, I thought "Oh, they found a good fix". But no! The values are "hardcoded". mtext(expression(list(mu[1]==0,mu[2]==0,sigma[11]==10,sigma[22]==10,sigma[12]==15,rho==0.5)), side=3) There should be substitute or bquote there to get the values of the mu's and sigma's that are defined above. I think I will register on your wiki and post a fix on that one, while I'm at it. > Romain > > Paul Johnson wrote: >> >> Greetings: >> >> I would like comments on this example and after fixing it up, I need >> help from someone who has access to insert this in R's help page for >> plotmath. >> >> I uploaded a drawing >> http://pj.freefaculty.org/R/Normal-2009.pdf >> >> that is created by the following code >> >> http://pj.freefaculty.org/R/Normal1_2009_plotmathExample.R >> >> This will be a good addition to the plotmath help page/example. >> >> Here's why: >> >> plotmath can be somewhat challenging. A few more concrete examples >> of ways in which users can combine values from R programs with >> plotmath expressions will be helpful to new users. In particular, >> combining several values from a program into a single expression is a >> difficult task that is not currently exemplified. I believe this is >> an instructive example because it combines the use of lines, arrows, >> text, and polygons. >> >> Here's the code: >> >> ###Set mu and sigma at your pleasure: >> mu <- 10.03 >> sigma <- 12.5786786 >> >> myx <- seq( mu - 3.5*sigma, mu+ 3.5*sigma, length.out=500) >> >> myDensity <- dnorm(myx,mean=mu,sd=sigma) >> >> >> # It is challenging to combine plotmath with values of mu and sigma in >> one expression. >> # Either bquote or substitute can be used. First use bquote: >> >> myTitle1 <- bquote (paste("x ~ Normal(", mu== .(round(mu,2)), ',', >> sigma== .(round(sigma,2)),")") ) >> >> ### Using substitute: >> ### myTitle1 <- substitute( "x ~ Normal" ~~ group( "(", list(mu==mu1, >> sigma^2==sigma2)#, ")") , list(mu1=round(mu,2), >> sigma2=round(sigma^2,2))) >> >> ### Or combine the two: >> ### t1 <- substitute ( mu == a , list (a = mu)) >> ### t2 <- substitute ( sigma == a, list(a = round(sigma,2))) >> ### myTitle1 <- bquote(paste("x ~ Normal(", .(t1),",", .(t2),")" ) ) >> >> >> ## To save the plot in a file, change FALSE to TRUE >> saveplot <- FALSE >> >> if (saveplot == TRUE){ >> pdf(file="Normal-2009.pdf", height=6.5, width=6.5, onefile=F, >> paper="special", pointsize=10) >> >> }else { >> dev.new(height=6.5, width=6.5,pointsize=9) >> } >> ### xpd needed to allow writing outside strict box of graph >> par(xpd=TRUE, ps=10) >> >> plot(myx, myDensity, type="l", xlab="x", ylab="Probability Density ", >> main=myTitle1, axes=FALSE) >> axis(2, pos= mu - 3.6*sigma) >> axis(1, pos=0) >> lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes >> >> # bquote creates an expression that text plotters can use >> t1 <- bquote( mu== .(mu)) >> >> # Find a value of myx that is "very close to" mu >> centerX <- max(which (myx <= mu)) >> # plot light vertical line under peak of density >> lines( c(mu, mu), c(0, myDensity[centerX]), lty= 14, lwd=.2) >> >> # label the mean in the bottom margin >> mtext(bquote( mu == .(mu)), 1, at=mu, line=-1) >> >> ### find position 20% "up" vertically, to use for arrow coordinate >> ss = 0.2 * max(myDensity) >> # Insert interval to represent width of one sigma >> arrows( x0=mu, y0= ss, x1=mu+sigma, y1=ss, code=3, angle=90, length=0.1) >> >> ### Write the value of sigma above that interval >> t2 <- bquote( sigma== .(round(sigma,2))) >> text( mu+0.5*sigma, 1.15*ss, t2) >> >> ### Create a formula for the Normal >> normalFormula <- expression (f(x) == frac (1, sigma* sqrt(2*pi)) * >> e^{~~ - ~~ frac(1,2)~~ bgroup("(", frac(x-mu,sigma),")")^2}) >> # Draw the Normal formula >> text ( mu + 0.5*sigma, max(myDensity)- 0.10 * max(myDensity), >> normalFormula, pos=4) >> >> ### Theory says we should have 2.5% of the area to the left of: -1.96 * >> sigma. >> ### Find the X coordinate of that "critical value" >> criticalValue <- mu -1.96 * sigma >> ### Then grab all myx values that are "to the left" of that critical >> value. >> specialX <- myx[myx <= criticalValue] >> >> ### mark the critical value in the graph >> text ( criticalValue, 0 , label= paste(round(criticalValue,2)), pos=1) >> ### Take sequence parallel to values of myx inside critical region >> specialY <- myDensity[myx < criticalValue] >> # Polygon makes a nice shaded illustration >> polygon(c(specialX[1], specialX, specialX[length(specialX )]), c(0, >> specialY, 0), density=c(-110),col="lightgray" ) >> >> shadedArea <- round(pnorm(mu - 1.96 * sigma, mean=mu, sd=sigma),4) >> >> >> ### I want to insert annotation about area on left side. >> >> al1 <- bquote(Prob(x <= .(round(criticalValue,2)))) >> al2 <- bquote(phantom(0) == F( .(round(criticalValue,2)) )) >> al3 <- bquote(phantom(0) == .(round(shadedArea,3))) >> >> ### Hard to position this text "just right" >> ### Have tried many ideas, this may be least bad. >> ### Get center position in shaded area >> medX <- median(specialX) >> ### density at that center point: >> denAtMedX <- myDensity[indexMed <- max(which(specialX < medX))] >> >> text(medX, denAtMedX+0.0055, labels=al1) >> text(medX, denAtMedX+0.004, labels=al2) >> text(medX, denAtMedX+0.0025, labels=al3) >> >> ### point from text toward shaded area >> arrows( x0=medX, y0=myDensity[indexMed]+0.002 ,x1= mu-2.5 *sigma, y1= >> 0.40*myDensity[length(specialX)] , length=0.1) >> >> >> >> ss <- 0.1 * max(myDensity) >> ### Mark interval from mu to mu-1.96*sigma >> arrows( x0=mu, y0= ss, x1=mu-1.96*sigma, y1=ss, code=3, angle=90, >> length=0.1) >> ### Put text above interval >> text( mu - 2.0*sigma,1.15*ss, >> bquote(paste(.(round(criticalValue,2)),phantom(1)==mu-1.96 * >> sigma,sep="")),pos=4 ) >> >> >> >> >> criticalValue <- mu +1.96 * sigma >> ### Then grab all myx values that are "to the left" of that critical >> value. >> specialX <- myx[myx >= criticalValue] >> >> ### mark the critical value in the graph >> text ( criticalValue, 0 , label= paste(round(criticalValue,2)), pos=1) >> ### Take sequence parallel to values of myx inside critical region >> specialY <- myDensity[myx > criticalValue] >> # Polygon makes a nice shaded illustration >> polygon(c(specialX[1], specialX, specialX[length(specialX )]), c(0, >> specialY, 0), density=c(-110),col="lightgray" ) >> >> shadedArea <- round(pnorm(mu + 1.96 * sigma, mean=mu, sd=sigma, >> lower.tail=F),4) >> >> >> ### Insert simpler comment on right side. >> >> al2 <- bquote(phantom(0) == 1 - F( .(round(criticalValue,2)) )) >> al3 <- bquote(phantom(0) == .(round(shadedArea,3))) >> >> ### Hard to position this text "just right" >> ### Have tried many ideas, this may be least bad. >> ### Get center position in shaded area >> medX <- median(specialX) >> ### density at that center point: >> denAtMedX <- myDensity[indexMed <- max(which(specialX < medX))] >> >> >> text(medX, denAtMedX+0.004, labels=al2) >> text(medX, denAtMedX+0.0025, labels=al3) >> >> ### point from text toward shaded area >> arrows( x0=medX, y0=myDensity[indexMed]+0.002 ,x1= mu+2.5 *sigma, y1= >> 0.40*myDensity[length(specialX)] , length=0.1) >> >> >> >> ss <- 0.05 * max(myDensity) >> ### Mark interval from mu to mu+1.96*sigma >> arrows( x0=mu, y0= ss, x1=mu+1.96*sigma, y1=ss, code=3, angle=90, >> length=0.1) >> ### Put text above interval >> text( mu + 1.96*sigma,1.15*ss, >> bquote(paste(.(round(criticalValue,2)),phantom(1)==mu+1.96 * >> sigma,sep="")),pos=2 ) >> >> >> if (saveplot == TRUE) dev.off() >> > > -- > Romain Francois > Independent R Consultant > +33(0) 6 28 91 30 30 > http://romainfrancois.blog.free.fr > > > > -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
