# Example of breaking out of a loop using return
# the 'break' keyword could also be used, but is not
# exactly equivalent
> f <- function() {}
> f <- function() {
+   for (i in 1:10)
+      return(i)
+ }

# Function to compute the length of a vector using
# a for loop.
> f()
[1] 1
> l <- function(v) {
+   len <- 0
+   for (i in v)
+     len <- len + 1
+   len
+ }
> l(1:10)
[1] 10

# duplicate function using a for loop
> dup <- function(x, n=2) {
+   res <- c()
+   for (i in 1:n)
+      res <- c(res, x)
+   res
+ }
> dup(1)
[1] 1 1
> dup(1, n=10)
 [1] 1 1 1 1 1 1 1 1 1 1

# a recursive implementation of the above dup function
> dup2 <- function(x, n=2) {
+   if (n > 1)
+     c(x, dup2(x, n-1))
+   else
+     x
+ }
> dup2(1)
[1] 1 1
> dup2(1, n=4)
[1] 1 1 1 1

# Make an x-axis for our plot.  We will evaluate the
# dnorm function at each of these x points.
> xs <- seq(-4, 4, by=0.01)
> xs
  [1] -4.00 -3.99 -3.98 -3.97 -3.96 -3.95 -3.94 -3.93 -3.92
 [10] -3.91 -3.90 -3.89 -3.88 -3.87 -3.86 -3.85 -3.84 -3.83
 [19] -3.82 -3.81 -3.80 -3.79 -3.78 -3.77 -3.76 -3.75 -3.74
...
[775]  3.74  3.75  3.76  3.77  3.78  3.79  3.80  3.81  3.82
[784]  3.83  3.84  3.85  3.86  3.87  3.88  3.89  3.90  3.91
[793]  3.92  3.93  3.94  3.95  3.96  3.97  3.98  3.99  4.00

# Plot a standard normal curve for the x values in xs
> plot(xs, dnorm(xs))

# Same plot, but with fewer points on the x axis
> xs <- seq(-4, 4, by=0.1)
> plot(xs, dnorm(xs))

# Change the plotting from points to lines
> plot(xs, dnorm(xs), type="l")

# The distribution of the null distribution in HW3, #3
> plot(xs, dnorm(xs, sd=10/sqrt(25)), type="l")

# Increase the boundaries of the x-axis
> xs <- seq(-10, 10, by=0.1)
> plot(xs, dnorm(xs, sd=10/sqrt(25)), type="l")
> plot(xs, dnorm(xs, mean=1.5, sd=10/sqrt(25)), type="l")
# Using the plot function again destroyed our old plot!

# Create the first curve with plot(), the second with lines()
> plot(xs, dnorm(xs, sd=10/sqrt(25)), type="l")
> lines(xs, dnorm(xs, mean=1.5, sd=10/sqrt(25)), type="l", lty="dashed")
# To differentiate H0's curve from H1, make H1's curve dashed

# Plot a vertical line at the one-sided critical value for alpha=0.05
> qnorm(0.95, sd=10/sqrt(25))
[1] 3.289707
> abline(v=qnorm(0.95, sd=10/sqrt(25)))
> crit <- qnorm(0.95, sd=10/sqrt(25))
# Plot the y-value of the normal curve at the critical value
> text(0, 0.07, format(dnorm(crit, sd=2), digits=3))

# Fill in the area under the curve
> polygon(c(crit, xs[xs >= crit]), c(0, dnorm(xs[xs >= crit], mean=1.5, sd=2)), col="grey")

# Same plot, but for the alternative hypothesis (mean=3)
> plot(xs, dnorm(xs, sd=10/sqrt(25)), type="l")
> lines(xs, dnorm(xs, mean=3, sd=10/sqrt(25)), type="l", lty="dashed")
> abline(v=crit)
> polygon(c(crit, xs[xs >= crit]), c(0, dnorm(xs[xs >= crit], mean=3, sd=2)), col="grey")


# Open up a new plotting window so that plot commands do not
# modify the previous plot.  'quartz' is Mac-specific.  Linux
# (and probably Windows) users should use x11().
# Here we plot the curve and power for the null hypothesis
# (the same plot that we showed previously, but we can now
# compare the two plots side by side).
> quartz()
> plot(xs, dnorm(xs, sd=10/sqrt(25)), type="l")
> lines(xs, dnorm(xs, mean=1.5, sd=10/sqrt(25)), type="l", lty="dashed")
> abline(v=qnorm(0.95, sd=10/sqrt(25)))
> polygon(c(crit, xs[xs >= crit]), c(0, dnorm(xs[xs >= crit], mean=1.5, sd=2)), col="grey")



# Plot a standard normal curve and show the function
# that defines area under the curve up to the point x.
> plot(xs, dnorm(xs), type='l')
> polygon(c(-10, xs[xs <= 1], 1), c(0, dnorm(xs[xs <= 1]), 0), col="grey")
> text(5, 0.2, quote(Phi(x) == frac(1, sqrt(2*pi)*sigma) * integral(e^frac(-(x - mu)^2, 2*sigma^2)*dx, -infinity, x)))
# The function works for all x.  The shaded area shows the
# function at x=1 (Phi(1)).

# Increase the size of the text using the cex parameter
> text(5, 0.2, quote(Phi(x) == frac(1, sqrt(2*pi)*sigma) * integral(e^frac(-(x - mu)^2, 2*sigma^2)*dx, -infinity, x)), cex=1.5)


> ?points
starting httpd help server ... done
# Copied directly from the points manual 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)
> demo(plotmath)


	demo(plotmath)
	---- ~~~~~~~~

Type  <Return>	 to start : 
...


# Create a matrix of normally distributed random values
> m <- matrix(rnorm(150), nrow=10)
> image(m)
> dim(m)
[1] 10 15

# Image plots m sideways.  Transpose and flip m to correct this.
> image(t(m[nrow(m):1,]))


# The layout function
> layout(matrix(1:4, ncol=2))  # Split the screen into 4 equal plots
> layout.show(4)               # Shows the boundaries and positions of each plot area
> layout(matrix(1:2, ncol=1))  # Split into two plots
> layout.show(2)
> layout(matrix(1:2, ncol=1), heights=matrix(c(3, 1), ncol=1))
> layout.show(2)
> layout(matrix(1:2, ncol=1), heights=matrix(c(3, 1), ncol=1))
> image(t(m[nrow(m):1,]))

# Building a color key
> range(m)                 # First take the range of values in the matrix
[1] -2.423600  2.895753
> r <- range(m)
# Make a sequence of values from r[1] to r[2] stepping by 0.01
> seq(r[1], r[2], by=0.01)
  [1] -2.423599516 -2.413599516
  [3] -2.403599516 -2.393599516
  [5] -2.383599516 -2.373599516
...
[529]  2.856400484  2.866400484
[531]  2.876400484  2.886400484
> cs <- seq(r[1], r[2], by=0.01)
> image(cs)
Error in image.default(cs) : argument must be matrix-like
# The image function requires matrices!
> image(as.matrix(cs))
Error in plot.new() : figure margins too large
# Due to the layout function, we no longer have enough room
# in our plotting area to contain the margins of the plot.

# Make the margins smaller.  This option will persist
# until the current plotting area is closed or it is
# manually changed again.
> par(mar=c(1, 1, 1, 1))

# Now our color key will not complain about not having
# enough room.
> image(as.matrix(cs))

# There wasn't enough room to see the x-axis of the color key!
# Make the bottom margin bigger.
> par(mar=c(2, 1, 1, 1))

# Start over and replot the heatmap and colorkey.  This
# time, yaxt="n" means the colorkey will not have a y-axis.
> image(t(m[nrow(m):1,]))
> image(as.matrix(cs), yaxt="n")

# Give y-values to each row instead of using the
# automatically assigned values.
> image(y=1:10, z=t(m[nrow(m):1,]), yaxt="n")
# Now we can plot labels for the rows at the correct
# positions.  las=2 means that the labels should always
# be oriented perpendicular to the axis.
> axis(side=2, at=1:10, labels=letters[1:10], las=2)
# Showing that our heatmap's y-axis has indeed changed
> text(0.1, 2, "test")

# Close the previous plot

# Finally, plot the heatmap and color key with axes using
# letters as row names (which will be perpindicular to
# the axis due to las=2), and without drawing tick marks.
> image(y=1:10, z=t(m[nrow(m):1,]), yaxt="n")
> axis(side=2, at=1:10, labels=letters[1:10], las=2, tick=FALSE)
> image(x=cs, z=as.matrix(cs), yaxt="n")


# Loading the gplots library for the colorpanel function
> library(gplots)
Loading required package: gtools
Loading required package: gdata

Attaching package: 'gdata'


	The following object(s) are masked from package:utils :

	 object.size 

Loading required package: caTools
Loading required package: bitops
Loading required package: grid

Attaching package: 'gplots'


	The following object(s) are masked from package:stats :

	 lowess 


# How one would install the gplots package if it weren't already installed.
# Requires an internet connection.
> install.packages("gplots")
trying URL 'http://cran.fhcrc.org/bin/macosx/leopard/contrib/2.10/gplots_2.8.0.tgz'
Content type 'application/x-gzip' length 333054 bytes (325 Kb)
opened URL
==================================================
downloaded 325 Kb


The downloaded packages are in
	/var/folders/T5/T5Z4ttKEHwSf84L9E0-kWE+++TM/-Tmp-//Rtmpe6yYwx/downloaded_packages

# The same heatmap as before but with a different color key
> image(y=1:10, z=t(m[nrow(m):1,]), yaxt="n", col=colorpanel(100, low="blue", mid="white", high="red"))


# An escape-velocity colored plot of the Mandelbrot set
> load("mandel.E.RData")  # Computed before class and saved using save()
# The matrix 'E' was contained in the RData file
> plot.mandel <- function(E) image(z=t(E), col=colorpanel(100, low="blue", mid="white", high="black"))
> plot.mandel(E)