inst/1-programming.R

In msme: Functions and Datasets for "Methods of Statistical Model Estimation"

### R code from vignette source '2-programming.rnw'

###################################################
### code chunk number 1: Setup
###################################################
options(repos="http://cran.r-project.org")

if(!require(Hmisc, quietly=TRUE)) install.packages("Hmisc")
if(!require(xtable, quietly=TRUE)) install.packages("xtable")

rm(list=ls())

options(width = 67)


###################################################
### code chunk number 2: 2-programming.rnw:104-105
###################################################
1 + 2


###################################################
### code chunk number 3: 2-programming.rnw:159-160
###################################################
wavelengths <- c(325.3, 375.6, 411.1)


###################################################
### code chunk number 4: 2-programming.rnw:169-170
###################################################
wavelengths


###################################################
### code chunk number 5: 2-programming.rnw:179-180
###################################################
class(wavelengths)


###################################################
### code chunk number 6: 2-programming.rnw:190-192
###################################################
sentence <- c("This", "is", "a", "character", "vector")
class(sentence)


###################################################
### code chunk number 7: 2-programming.rnw:203-204
###################################################
wavelengths / 1000


###################################################
### code chunk number 8: 2-programming.rnw:209-210
###################################################
mean(wavelengths)


###################################################
### code chunk number 9: 2-programming.rnw:214-215
###################################################
length(wavelengths)


###################################################
### code chunk number 10: 2-programming.rnw:227-229
###################################################
a_factor <- factor(c("A","A","B","B","B","C"))
a_factor


###################################################
### code chunk number 11: 2-programming.rnw:235-237
###################################################
a_factor[6] <- "Z"
a_factor


###################################################
### code chunk number 12: 2-programming.rnw:248-249
###################################################
levels(a_factor)


###################################################
### code chunk number 13: 2-programming.rnw:252-253
###################################################
levels(a_factor)[2] <- "Bee"


###################################################
### code chunk number 14: 2-programming.rnw:257-258
###################################################
a_factor


###################################################
### code chunk number 15: 2-programming.rnw:271-272
###################################################
factor(1) + factor(2)


###################################################
### code chunk number 16: 2-programming.rnw:281-282
###################################################
TRUE + TRUE


###################################################
### code chunk number 17: 2-programming.rnw:289-290
###################################################
1:4 < 3


###################################################
### code chunk number 18: 2-programming.rnw:298-299
###################################################
TRUE | TRUE


###################################################
### code chunk number 19: 2-programming.rnw:316-318
###################################################
missing.bits <- c(1, NA, 2)
mean(missing.bits)


###################################################
### code chunk number 20: 2-programming.rnw:324-325
###################################################
mean(missing.bits, na.rm = TRUE)


###################################################
### code chunk number 21: 2-programming.rnw:329-330
###################################################
length(missing.bits)


###################################################
### code chunk number 22: 2-programming.rnw:348-349
###################################################
sentence[2]


###################################################
### code chunk number 23: 2-programming.rnw:353-354
###################################################
sentence[1:3]


###################################################
### code chunk number 24: 2-programming.rnw:360-362
###################################################
sentence[-4]
sentence[-(2:4)]


###################################################
### code chunk number 25: 2-programming.rnw:366-368
###################################################
wavelengths > 400
wavelengths[wavelengths > 400]


###################################################
### code chunk number 26: 2-programming.rnw:391-392
###################################################
my.list <- list(number = 1, text = "alphanumeric")


###################################################
### code chunk number 27: 2-programming.rnw:399-403
###################################################
my.list[[1]]
my.list$text
my.list[[1]] <- 2
my.list


###################################################
### code chunk number 28: 2-programming.rnw:408-409
###################################################
capacious.empty.list <- vector(1000, mode = "list")


###################################################
### code chunk number 29: 2-programming.rnw:419-420
###################################################
sapply(my.list, class)


###################################################
### code chunk number 30: 2-programming.rnw:446-449
###################################################
example <- data.frame(var.a = 1:3,
                      var.b = c("a","b","c"))
str(example)


###################################################
### code chunk number 31: 2-programming.rnw:457-459
###################################################
data.frame(var.a = 1,
           var.b = c("a","b","c"))


###################################################
### code chunk number 32: 2-programming.rnw:474-475
###################################################
example[2, 1:2]


###################################################
### code chunk number 33: 2-programming.rnw:479-480
###################################################
example[2,]


###################################################
### code chunk number 34: subset
###################################################
subset(example, subset = var.a > 1, select = "var.b")


###################################################
### code chunk number 35: 2-programming.rnw:519-525
###################################################
example.ok <- function (a, b = 1) {
  return(a + b)
}

example.ok(2,2)
example.ok(2)


###################################################
### code chunk number 36: 2-programming.rnw:532-533
###################################################
formals(example.ok)     


###################################################
### code chunk number 37: 2-programming.rnw:536-537
###################################################
body(example.ok)        


###################################################
### code chunk number 38: 2-programming.rnw:540-541
###################################################
environment(example.ok) 


###################################################
### code chunk number 39: 2-programming.rnw:565-567
###################################################
example <- function(a = 1) a
example()


###################################################
### code chunk number 40: 2-programming.rnw:574-577
###################################################
example <- function(a, b) a
example(a = 1)
example(a = 1, b = log(-1))


###################################################
### code chunk number 41: 2-programming.rnw:589-590
###################################################
example <- function(a) a + b


###################################################
### code chunk number 42: 2-programming.rnw:603-605
###################################################
b <- 1
example(a = 1)


###################################################
### code chunk number 43: 2-programming.rnw:616-621
###################################################
example <- function(a, b) a + b

contains <- function(a, ...) example(a, ...)

contains(a = 1, b = 2)


###################################################
### code chunk number 44: 2-programming.rnw:629-631
###################################################
example(1, 2)
example(1, b = 2)


###################################################
### code chunk number 45: 2-programming.rnw:654-659
###################################################
a <- 1
a
risky <- function(x) a <<- x
risky(2)
a


###################################################
### code chunk number 46: 2-programming.rnw:667-671
###################################################
example.ok <- function (a, b = 1) {
  a + b
}
example.ok(2)


###################################################
### code chunk number 47: 2-programming.rnw:681-685
###################################################
example.ok <- function (a, b = 1) {
  return(list(a = a, b = b, sum = a + b))
}
example.ok(2)


###################################################
### code chunk number 48: 2-programming.rnw:690-692
###################################################
example.ok(2)$sum
example.ok(2)[[3]]


###################################################
### code chunk number 49: 2-programming.rnw:713-716
###################################################
x <- 100
scope.fn <- function() x / 10
scope.fn()


###################################################
### code chunk number 50: 2-programming.rnw:738-739
###################################################
variance.binomial <- function(mu) mu * (1 - mu / m) 


###################################################
### code chunk number 51: 2-programming.rnw:744-748
###################################################
irls <- function(mu) {
  m <- 100
  variance.binomial(mu)
}


###################################################
### code chunk number 52: 2-programming.rnw:767-768
###################################################
environment(variance.binomial)


###################################################
### code chunk number 53: 2-programming.rnw:775-781
###################################################
irls <- function(mu) {
  m <- 100
  variance.binomial(mu, m)
}
variance.binomial <- function(mu, m) mu * (1 - mu / m) 
irls(0.5)


###################################################
### code chunk number 54: 2-programming.rnw:786-792
###################################################
irls <- function(mu) {
   variance.binomial <- function(mu) mu * (1 - mu / m)
   m <- 100
   variance.binomial(mu)
   }
irls(0.5)


###################################################
### code chunk number 55: 2-programming.rnw:829-831
###################################################
mat1 <- matrix(1:4, nrow=2)
mat1


###################################################
### code chunk number 56: 2-programming.rnw:838-840
###################################################
attributes(mat1)
dim(mat1)


###################################################
### code chunk number 57: 2-programming.rnw:852-854
###################################################
attr(wavelengths, "units") <- "micrometres"
attributes(wavelengths)


###################################################
### code chunk number 58: 2-programming.rnw:862-863
###################################################
sum(mat1)


###################################################
### code chunk number 59: 2-programming.rnw:871-872
###################################################
colSums(mat1)


###################################################
### code chunk number 60: 2-programming.rnw:881-884
###################################################
(mat2 <- diag(2))
mat1 + mat2
mat1 * mat2


###################################################
### code chunk number 61: 2-programming.rnw:892-893
###################################################
mat1 %*% mat2


###################################################
### code chunk number 62: 2-programming.rnw:902-903
###################################################
tcrossprod(mat1, mat2)


###################################################
### code chunk number 63: 2-programming.rnw:909-910
###################################################
solve(mat1)


###################################################
### code chunk number 64: 2-programming.rnw:913-914
###################################################
mat1 %*% solve(mat1)


###################################################
### code chunk number 65: 2-programming.rnw:967-968
###################################################
dnorm(0, mean = 0, sd = 1)


###################################################
### code chunk number 66: 2-programming.rnw:972-973
###################################################
dnorm(0)


###################################################
### code chunk number 67: 2-programming.rnw:988-989
###################################################
pnorm(0.5)


###################################################
### code chunk number 68: 2-programming.rnw:994-995
###################################################
qnorm(0.975, mean = 2, sd = 2)


###################################################
### code chunk number 69: 2-programming.rnw:1000-1001
###################################################
rnorm(5, mean = 2, sd = 2)


###################################################
### code chunk number 70: 2-programming.rnw:1023-1026
###################################################
dwatson <- function(x, theta) {
  (1 + theta) / theta / (1 + x / theta)^2
}


###################################################
### code chunk number 71: 2-programming.rnw:1049-1053
###################################################
pwatson <- function(q, theta) {
  integrate(function(x) dwatson(x, theta), 
            lower = 0, upper = q)$value
  }


###################################################
### code chunk number 72: 2-programming.rnw:1064-1068
###################################################
qwatson <- function(p, theta) {
  uniroot(function(x) pwatson(x, theta) - p, 
          lower = .Machine$double.eps, upper = 1)$root
}


###################################################
### code chunk number 73: 2-programming.rnw:1086-1087
###################################################
pwatson(1, theta = 1)


###################################################
### code chunk number 74: 2-programming.rnw:1099-1102
###################################################
rwatson <- function(n, theta) {
  sapply(runif(n), function(x) qwatson(x, theta))
}


###################################################
### code chunk number 75: 2-programming.rnw:1106-1107
###################################################
rwatson(5, 1)


###################################################
### code chunk number 76: compare
###################################################
system.time(rnorm(1000))
system.time(rwatson(1000,1))


###################################################
### code chunk number 77: 2-programming.rnw:1125-1132
###################################################
qwatson1 <- function(p, theta) {
  (1 / ( 1 - p / (1 + theta) ) - 1) * theta
}
rwatson1 <- function(n, theta) {
  qwatson1(runif(n), theta)
}
system.time(rwatson1(1000,1))


###################################################
### code chunk number 78: 2-programming.rnw:1224-1225
###################################################
for (i in 1:2) cat(i, "\n")


###################################################
### code chunk number 79: 2-programming.rnw:1230-1231
###################################################
for (i in c("a","b")) cat(i, "\n")


###################################################
### code chunk number 80: 2-programming.rnw:1236-1237
###################################################
i


###################################################
### code chunk number 81: 2-programming.rnw:1246-1248
###################################################
seq(2)
seq(NULL)


###################################################
### code chunk number 82: 2-programming.rnw:1259-1262
###################################################
output <- vector(length = 10, mode = "list")
for (i in 1:10) output[[i]] <- paste("Element", i)
output[[1]]


###################################################
### code chunk number 83: 2-programming.rnw:1268-1270
###################################################
output <- lapply(1:10, function(i) paste("Element", i))
output[[1]]


###################################################
### code chunk number 84: 2-programming.rnw:1286-1291
###################################################
i <- 1
while (i < 2) {
  cat(i, "\n")
  i <- i + 1
}


###################################################
### code chunk number 85: fig-test-traj
###################################################
test.fun <- function(x) 2*x^2 - log(3*x)
curve(test.fun, from = 0, to = 3)


###################################################
### code chunk number 86: test-traj
###################################################
par(las=1, mar=c(4,4,2,1))
test.fun <- function(x) 2*x^2 - log(3*x)
curve(test.fun, from = 0, to = 3)


###################################################
### code chunk number 87: 2-programming.rnw:1546-1547
###################################################
options(error = recover)


###################################################
### code chunk number 88: 2-programming.rnw:1590-1591
###################################################
length(apropos("^print."))


###################################################
### code chunk number 89: 2-programming.rnw:1596-1597
###################################################
sessionInfo()[[1]]$version.string


###################################################
### code chunk number 90: 2-programming.rnw:1629-1630
###################################################
print


###################################################
### code chunk number 91: 2-programming.rnw:1655-1657
###################################################
ordinal <- ordered(c(1,2,3))
class(ordinal)


###################################################
### code chunk number 92: 2-programming.rnw:1668-1669
###################################################
ordinal


###################################################
### code chunk number 93: 2-programming.rnw:1678-1680
###################################################
class(ordinal) <- "ordered"
ordinal


###################################################
### code chunk number 94: 2-programming.rnw:1690-1691
###################################################
methods(class = "ordered")


###################################################
### code chunk number 95: 2-programming.rnw:1699-1700
###################################################
methods(nobs)


###################################################
### code chunk number 96: 2-programming.rnw:1733-1735
###################################################
item <- 1:2
class(item) <- c("thing2","thing1")


###################################################
### code chunk number 97: 2-programming.rnw:1739-1742
###################################################
print.thing1 <- function(x, ...) {
  cat("inherits from Thing 1.\n")
}


###################################################
### code chunk number 98: 2-programming.rnw:1755-1759
###################################################
print.thing2 <- function(x, ...) {
  cat("Thing 2 ")
  NextMethod()
}


###################################################
### code chunk number 99: 2-programming.rnw:1763-1764
###################################################
item


###################################################
### code chunk number 100: 2-programming.rnw:1769-1773
###################################################
print.thing3 <- function(x, ...) {
  cat("Thing 3 also ")
  NextMethod()
}


###################################################
### code chunk number 101: 2-programming.rnw:1777-1780
###################################################
another.item <- 1:3
class(another.item) <- c("thing3","thing1")
another.item


###################################################
### code chunk number 102: 2-programming.rnw:1962-1963
###################################################
M <- matrix(c(3,4,6,8), nrow=1)


###################################################
### code chunk number 103: 2-programming.rnw:1966-1967
###################################################
P <- matrix(c(3,4,6,8,4,8,4,7,2,2,5,4,4,7,5,2), ncol=4)


###################################################
### code chunk number 104: 2-programming.rnw:1997-2040 (eval = FALSE)
###################################################
## available <- "I'm in the enclosure!\n"
## 
## test.ok <- function (a, b = 1) {
##   cat(available)
##   a + b
## }
## 
## test.ok(2)
## 
## test.ok <- function (a, ...) {
##   a + b
## }
## 
## test.ok(2)
## 
## test.still.ok <- function (a, b = 1) {
##   cat(still.available)
##   a + b
## }
## 
## still.available <- "I'm still in the enclosure!\n"
## 
## test.still.ok(2)
## 
## test.not.ok <- function(a, ...) {
##   available <- "I'm in the enclosure!\n"
##   test.ok <- function (a, ...) {
##     cat(available)
##     a + b
##   }
##   print(test.ok(2))  
##   }
## 
## test.not.ok <- function(a, ...) {
##   test.ok <- function (a, ...) {
##     a + b
##   }
##   print(test.ok(1, b = 2))  
##   print(test.ok(1, ...))  
##   }
## 
## test.not.ok(1)
##