R devel - Jun 2009 - Would like to add this to example for plotmath. Can you help?

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,
y10.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,
y10.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()


-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas

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