## Settings for RMarkdown http://yihui.name/knitr/options#chunk_options
opts_chunk$set(comment = "", warning = FALSE, message = FALSE, tidy = FALSE,
echo = TRUE, fig.width = 10, fig.height = 8)
options(width = 116, scipen = 10)
setwd("~/statistics/bio509/")
Task 1
## Define a recursive function
legendre <- function(n, x) {
## Constant part
if (n == 0) return(1)
if (n == 1) return(x)
## Recursion part
((2 * n - 1) * x * Recall(n - 1, x) - (n - 1) * Recall(n - 2, x)) / n
}
## Vectorize
legendre.vec <- Vectorize(legendre, c("n","x"))
## Result
mat.legendre <- outer(X = 0:6, Y = 0:6, FUN = function(X, Y) legendre.vec(x = X, n = Y))
colnames(mat.legendre) <- paste("n=", 0:6, sep = "")
rownames(mat.legendre) <- paste("x=", 0:6, sep = "")
mat.legendre
n=0 n=1 n=2 n=3 n=4 n=5 n=6
x=0 1 0 -0.5 0 0.375 0.0 -0.3125
x=1 1 1 1.0 1 1.000 1.0 1.0000
x=2 1 2 5.5 17 55.375 185.8 634.9375
x=3 1 3 13.0 63 321.000 1683.0 8989.0000
x=4 1 4 23.5 154 1060.375 7511.5 54200.6875
x=5 1 5 37.0 305 2641.000 23525.0 213445.0000
x=6 1 6 53.5 531 5535.375 59357.2 648316.9375
## Confirmation
n0 <- function(x) 1
n1 <- function(x) x
n2 <- function(x) 1/2 * (3 * x^2 - 1)
n3 <- function(x) 1/2 * (5 * x^3 - 3 * x)
n4 <- function(x) 1/8 * (35 * x^4 - 30 * x^2 + 3)
n5 <- function(x) 1/8 * (63 * x^5 - 70 * x^3 +15 * x)
n6 <- function(x) 1/16 * (231 * x^6 - 315 * x^4 + 105 * x^2 - 5)
data.frame(n0(0:6), n1(0:6), n2(0:6), n3(0:6), n4(0:6), n5(0:6), n6(0:6))
n0.0.6. n1.0.6. n2.0.6. n3.0.6. n4.0.6. n5.0.6. n6.0.6.
1 1 0 -0.5 0 0.375 0.0 -0.3125
2 1 1 1.0 1 1.000 1.0 1.0000
3 1 2 5.5 17 55.375 185.8 634.9375
4 1 3 13.0 63 321.000 1683.0 8989.0000
5 1 4 23.5 154 1060.375 7511.5 54200.6875
6 1 5 37.0 305 2641.000 23525.0 213445.0000
7 1 6 53.5 531 5535.375 59357.2 648316.9375
Task 2
plot(0, 0, type = "n", xlab = "x", ylab = "")
sequence <- seq(-10, 10, 0.001)
for (i in 0:6) {
lines(x = sequence, y = legendre.vec(n = i, x = sequence), col = i + 1)
}
legend("bottomright", legend = 0:6, col = 0:6 + 1, lty = 1, bty = "n")
## library(ggplot2)
## ggplot(data.frame(x = 0), aes(x)) + xlim(-5,5) + ylim(-5,5) +
## stat_function(fun = function(x) legendre.vec(n = 0, x = x), aes(color = "zero")) +
## stat_function(fun = function(x) legendre.vec(n = 1, x = x)) +
## stat_function(fun = function(x) legendre.vec(n = 2, x = x)) +
## stat_function(fun = function(x) legendre.vec(n = 3, x = x)) +
## stat_function(fun = function(x) legendre.vec(n = 4, x = x)) +
## stat_function(fun = function(x) legendre.vec(n = 5, x = x)) +
## stat_function(fun = function(x) legendre.vec(n = 6, x = x))
Reference http://rosettacode.org/wiki/Towers_of_Hanoi#R
hanoimove <- function(ndisks, from, to, via) {
if ( ndisks == 1 )
cat("move disk", "from peg", from, "to peg", to, "\n")
else {
Recall(ndisks-1, from, via, to)
Recall(1, from, to, via)
Recall(ndisks-1, via, to, from)
}
}
hanoimove(4,"A","B","C")
move disk from peg A to peg C
move disk from peg A to peg B
move disk from peg C to peg B
move disk from peg A to peg C
move disk from peg B to peg A
move disk from peg B to peg C
move disk from peg A to peg C
move disk from peg A to peg B
move disk from peg C to peg B
move disk from peg C to peg A
move disk from peg B to peg A
move disk from peg C to peg B
move disk from peg A to peg C
move disk from peg A to peg B
move disk from peg C to peg B