From: Paul Johnson <pauljohn32_at_gmail.com>

Date: Wed, 03 Jun 2009 16:52:07 -0500

Date: Wed, 03 Jun 2009 16:52:07 -0500

2009/6/3 Romain François <francoisromain_at_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_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-develReceived on Wed 03 Jun 2009 - 21:55:05 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Wed 03 Jun 2009 - 22:34:56 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*