# Visual overview of power (HW3, #3)
xs <- seq(-4, 4, by=0.01)
plot(xs, dnorm(xs, sd=2))            # Plots points!  Can we get a nice line?
plot(xs, dnorm(xs, sd=2), type="l")  # Type = "lines"

# Make room for the alternative distribution
xs <- seq(-10, 10, by=0.01)
plot(xs, dnorm(xs, sd=2), type="l")  # How do we draw another line on the same plot?
lines(xs, dnorm(xs, mean=1.5, sd=2), lty="dashed")

# Draw a vertical line at the critical value for a one-sided test, alpha=0.05
crit <- qnorm(0.95, sd=2)
abline(v=crit, lty="dashed", col="red")

# What's the power?  Area under the curve right of the critical value
polygon(c(crit, xs[xs >= crit]), c(0, dnorm(xs[xs >= crit], mean=1.5, sd=2)), col="grey")

# The power is..
1 - pnorm(crit, mean=1.5, sd=2)
text(4.5, 0.025, "0.185")            # Plot the number on the graph

# What happened when we shifted the mean?
xs <- seq(-10, 10, by=0.01)
plot(xs, dnorm(xs, sd=2), type="l")  # How do we draw another line on the same plot?
lines(xs, dnorm(xs, mean=3, sd=2), lty="dashed")
crit <- qnorm(0.95, sd=2)            # Critical value is still the same
abline(v=crit, lty="dashed", col="red")
# Added (10, 0) to the end to correct the polygon.
polygon(c(crit, xs[xs >= crit], 10), c(0, dnorm(xs[xs >= crit], mean=3, sd=2), 0), col="grey")
1 - pnorm(crit, mean=3, sd=2)
text(5, 0.035, "0.442")  

# Look at the two plots side by side


# Math annotation
x.axis <- seq(-4, 4, by=0.01)
plot(x.axis, dnorm(x.axis), type="l")
text(-2.5, 0.3, quote(Phi(x) == frac(1, sqrt(2*pi)*sigma) * integral(e^frac(-(x-mu)^2, 2*sigma^2)*dx, -infinity, x)))
lines(c(-2, -1), c(0.275, dnorm(-1)))

# More math annotation
?points   # To find the TestChars function
# From the 'points' documentation page
TestChars <- function(sign=1, font=1, ...)
{
   if(font == 5) { sign <- 1; r <- c(32:126, 160:254)
   } else if (l10n_info()$MBCS) r <- 32:126 else r <- 32:255
   if (sign == -1) r <- c(32:126, 160:255)
   par(pty="s")
   plot(c(-1,16), c(-1,16), type="n", xlab="", ylab="",
        xaxs="i", yaxs="i")
   grid(17, 17, lty=1)
   for(i in r) try(points(i%%16, i%/%16, pch=sign*i, font=font,...))
}
TestChars(font=5)

# A demo of many different notations
?plotmath
demo(plotmath)





# Multiple plots in a single figure: e.g., homework 3
power <- function(alpha=0.05, n=25, mu=1.5, sigma=10) {
    crit <- qnorm(1-alpha, mean=0, sd=sigma/sqrt(n))
    1 - pnorm(crit, mean=mu, sd=sigma/sqrt(n))
}
alpha <- seq(0, 1, by=0.01)
ns <- 1:125
sds <- seq(0, 50, by=0.1)
mus <- seq(0, 20, by=0.1)

layout(matrix(1:4, ncol=2))
par(cex=0.60)
plot(ns, power(n=ns), type="l", main="Power vs. sample size", xlab="Sample size", ylab="Power")
plot(alpha, power(alpha=alpha), type="l", main="Power vs. significance level", xlab="Significance level", ylab="Power")
plot(sds, power(sigma=sds), type="l", main="Power vs. standard deviation", xlab="Standard deviation", ylab="Power")
plot(mus, power(mu=mus), type="l", main="Power vs. mean difference", xlab="Difference between sample means", ylab="Power")





# Heatmaps and lots of common graphics options
# Package installation
install.packages("gplots")

# Heatmap on a random matrix
m <- matrix(rnorm(150), nrow=15)
image(z=m)

# Quirk: R plots the transpose and horizontal reflection of the matrix
image(z=t(m)[nrow(m):1,])

# Plotting a color key
layout(matrix(1:2, ncol=1), heights=matrix(c(2,1), ncol=1))
image(z=t(m[nrow(m):1,]))
x <- seq(r[1], r[2], by=0.01)
image(z=as.matrix(x))

# Big ugly space between the plots!  Margins
layout(matrix(1:2, ncol=1), heights=matrix(c(4,1), ncol=1))
image(z=t(m[nrow(m):1,]))
image(z=as.matrix(x)))
# Margin error! Can't just increase the relative sizes...
# Need to reduce the margins.
par(mar=c(5, 4, 0, 2))
image(z=as.matrix(seq(r[1], r[2], by=0.01)))

# y-axis tick marks on the color key don't make sense
# x-axis labels don't make sense
layout(matrix(1:2, ncol=1), heights=matrix(c(4,1), ncol=1))
image(z=t(m[nrow(m):1,]))
par(mar=c(5, 4, 0, 2))
image(x=x, z=as.matrix(x), yaxt="n")

# axis marks in general don't make sense.  Draw our own
layout(matrix(1:2, ncol=1), heights=matrix(c(4,1), ncol=1))
image(y=1:nrow(m), z=t(m[nrow(m):1,]), yaxt="n")
axis(side=2, at=1:nrow(m), labels=letters[1:nrow(m)])
# But... the text is rotated.  Make it perpendicular to the axis
layout(matrix(1:2, ncol=1), heights=matrix(c(4,1), ncol=1))
image(y=1:nrow(m), z=t(m[nrow(m):1,]), yaxt="n")
axis(side=2, at=nrow(m):1, labels=letters[1:nrow(m)], las=2)

# We might also remove the tick marks..
layout(matrix(1:2, ncol=1), heights=matrix(c(4,1), ncol=1))
image(y=1:nrow(m), z=t(m[nrow(m):1,]), yaxt="n")
axis(side=2, at=nrow(m):1, labels=letters[1:nrow(m)], las=2, tick=FALSE)
par(mar=c(5, 4, 0, 2))
image(x=x, z=as.matrix(x), yaxt="n")




# Mandelbrot set (general image plotting with image())
# Cairo library
library(Cairo)
compute.mandel <- function(g=1000, iter=25) {
        z <- matrix(0, nrow=g+1, ncol=1.5*g+1)
        res <- 2*(col(z)-1) / g - 2
        ims <- -(2 * (row(z) - 1) / g - 1)
	# P should be computed using outer()
        P <- matrix(complex(real=res, imag=ims), nrow=g+1, ncol=1.5*g+1)
        E <- matrix(Inf, nrow=g+1, ncol=1.5*g + 1)

        for (i in 1:iter) {
                z <- z^2 + P
                E.cur <- E
                E.cur[as.numeric(abs(z)) > 2] <- i
                E <- pmin(E, E.cur)
        }
        E[E == Inf] <- iter + 1
        E
}

plot.mandel <- function(E)
        image(z=t(E), col=colorpanel(100, low="blue", mid="white", high="black"))

# Not doing this: E <- compute.mandel()
load("mandel.E.RData")   # Load the escape metric from a previous R session
plot.mandel(E)
# Displaying the plot is somewhat quirky.  We could instead
# plot it to a PNG using the Cairo library.
png("mandel.png", width=1500, height=1000)  # Dimensions in pixels
plot.mandel(E)
dev.off()                                   # Complete the PNG