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