[R] Plotmath with values?

[Paul Johnson]

On Wed, Dec 31, 2008 at 10:47 AM, Mike Prager <mike.prager at noaa.gov> wrote:
> I hope to use the plotmath facility to print titles that mix
> math and values of R variables.
>
> The help for "plotmath" has an example, which after repeated
> reading, I find baffling. Likewise, I have read the help file
> for "substitute" (wqhich seems to be needed) without ever
> understanding what it does, other than being used in some magic
> incantations.
>
> I would like to do something like this:
>
>     dev.new()
>     aa <- round(pi,2)
>     plot(1:3, 1:3, main = ~ a == aa)
>
> and have the main title read "a = 3.14"
>
> but of course it reads "a = aa".
>

I agree this is a very difficult part of using R.  I asked this exact
same kind of question last year.  If you go to an R-help email archive
for April 2, 2008

"Trouble combining plotmath, bquote, expressions"

you will find some discussion in the list.

For a while, I got pretty excited and started working on better
example code to put into the R distribution, but lost initiative when
I saw the magnitude of the problem.  This example code will show
several usages of bquote and an alternative substitute approach.  The
result seems not to look so beautiful as I recall, but isn't that
always the way it goes :).


### Filename: Normal1_2008.R
### Paul Johnson March 31, 2008
### This code should be available somewhere in
http://pj.freefaculty.org.  If it is not
### email me <pauljohn at ku.edu>



mu <- 10.034

sigma <- 12.5786786

myx <- seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)

myDensity <- dnorm(myx,mean=mu,sd=sigma)


# Here's one way to retrieve the values of mu and sigma and insert
them into an expression
t1t2 <- bquote (paste("Normal(", mu== .(round(mu,2)), ',', sigma==
.(round(sigma,2)),")") )

plot(myx, myDensity, type="l", xlab="x", ylab="Probability Density ", main=t1t2)


t1t2 <- bquote (paste("Normal (   ", mu== .(round(mu,2)), '  ,  ',
sigma^2== .(round(sigma^2,2))," )") )
### Note spaces manually inserted above are needed, otherwise plotmath
overlaps "l" of normal with first parenthesis
plot(myx, myDensity, type="l", xlab="x", ylab="Probability Density ", main=t1t2)




### Had difficulty finding syntax for substitute to combine symbols mu
and sigma with values.
### Following is best I can figure, no simpler or obvious than previous method.
##t1 <- substitute ( mu == a ,  list (a = mu))
##t2 <- substitute ( sigma == a, list(a = sigma))
##t1t2 <- bquote(paste("Normal(", .(t1),",", .(t2),")" ) )



t1t2 <-  substitute( "Normal" ~~ group( "(", list(mu==mu1,
sigma^2==sigma2), ")") ,  list(mu1=round(mu,2),
sigma2=round(sigma^2,2)))


plot(myx, myDensity, type="l", xlab="x", ylab="Probability Density ", main=t1t2)



plot(myx, myDensity, type="l", xlab="x", ylab="Probability Density ",
main=t1t2, axes=F)
axis(2, pos= mu - 3.6*sigma)
axis(1, pos=0)


# bquote creates an expression that text plotters can use
t1 <-  bquote( mu== .(mu))

text( mu, 0.00, t1, pos=3)

ss = 0.2 * max(myDensity)

arrows( x0=mu, y0= ss, x1=mu+sigma, y1=ss, code=3, angle=90, length=0.1)

t2 <-  bquote( sigma== .(round(sigma,2)))

text( mu+0.5*sigma, 1.15*ss, t2)


normalFormula <- expression (f(x) == frac (1, sqrt(2*pi)*sigma) *
e^{~~ - ~~ frac(1,2)~~ bgroup("(", frac(x-mu,sigma),")")^2})

text ( mu + 0.5*sigma, max(myDensity)- 0.10 * max(myDensity),
normalFormula, pos=4)

### Theory says we should have about 2.5% of the area to the left of:
-1.96 * std.dev

criticalValue <- mu -1.96 * sigma
specialX <-  myx[myx <= criticalValue]

### mark the critical value in the graph

text ( criticalValue, 0 , label= paste(round(criticalValue,2)), pos=1)

specialY <- myDensity[myx < criticalValue]

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)

### Hard to position this text "just right"
al <- paste( "Prob(", "x <=" , round(criticalValue,2),")\n=","F(",
round(criticalValue,2) ,")\n=", round(shadedArea,3),sep="")
text(  criticalValue- sigma, myDensity[length(specialX)], labels=al, pos=3)

ss <- 0.1 * max(myDensity)
arrows( x0=mu, y0= ss, x1=mu-1.96*sigma, y1=ss, code=3, angle=90, length=0.1)

text( mu - 2.0*sigma, 1.15*ss,
bquote(paste(.(round(criticalValue,2)),phantom(1) == mu - 1.96," ",
sigma,sep=" ")),pos=4 )





> >From a user's point of view -- one who has never written a
> parser nor taken a course in compilers -- what is needed is the
> nonexistent function "value" usable in plotmath expressions to
> produce the value of its argument, as
>
>     plot(1:3, 1:3, main = ~ a == value(aa))
>
> How can this be done?
>
> THANKS!
>
> --
> Mike Prager, NOAA, Beaufort, NC
> * Opinions expressed are personal and not represented otherwise.
> * Any use of tradenames does not constitute a NOAA endorsement.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



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