## Chapter 7
## November 17, 2007
## Alan T. Arnholt
#
options(width=65)
library(PASWR)
#
# Example 7.4
m <- 50
coeff <- array(0, m)
for (n in 2:m)
{ coeff[n] <- (sqrt(2/(n-1))*gamma(n/2))/gamma((n-1)/2)}
plot(coeff[-1],type="l",xlab="n", ylab="coef",lwd=2)
abline(h=1,lty=2)
coeff[20]
# Or
curve(sqrt(2/(x-1))*gamma(x/2)/gamma((x-1)/2),2,50)
abline(h=1,lty=2)
#
# Robust Estimators Section:7.2.5
# Example 7.9
stem1 <- c(1.7, 2.8, 3.2, 3.4, 5.3, 5.9, 6.2, 7.2, 8.3, 9.3)
stem2 <- c(1.7, 2.8, 3.2, 3.4, 5.3, 5.9, 6.2, 7.2, 83, 9.3)
c(mean(stem1), sqrt(var(stem1)))
c(mean(stem2), sqrt(var(stem2)))
c(median(stem1), mad(stem1, constant = 1))
c(median(stem2), mad(stem2, constant = 1))
median(abs(stem1 - median(stem1)))
median(abs(stem2 - median(stem2)))
#
# Example 7.16 Oriental Cockroaches
attach(Roacheggs)
str(Roacheggs)      # Note: str(object) only works in R
mean(eggs)
p <- seq(0.1, 0.9, 0.001)
negloglike <- function(p){-(sum(eggs)*log(p) + sum(1-eggs)*log(1-p))}
nlm(negloglike, 0.2)
# nlmin(negloglike, 0.2)  # for S-PLUS
# Figure 7.3 (almost)
par(pty = "s")
p <- seq(0.1, 0.9, 0.001)
plot(p, - negloglike(p), type = "n", ylab = "L")
abline(v = mean(eggs), col = 13, lwd = 3)
lines(p, - negloglike(p), col = 6, lwd = 3)
#
# Optimize
loglike <- function(p){(sum(eggs)*log(p) + sum(1-eggs) * log(1-p))}
optimize(f=loglike,interval=c(0,1),maximum=TRUE)
detach(Roacheggs)
#
# Example 7.17
set.seed(23)
sum(rbinom(1000, 3, 0.5))/(1000 * 3)
#
set.seed(23)
mean(rbinom(1000, 3, 0.5))/3
#
# Example 7.18
set.seed(99)
mean(rpois(1000, 5))
#
# Example 7.20
# part c.
# generic Figure 7.4
par(pty = "s")
p <- seq(0.01, 0.99, 0.001)
loglike <- function(p)
{40 * log(2) - 220 * log(p) - 158 * log(1 - p)}
plot(p, - loglike(p), type = "n")
lines(p, - loglike(p), col = 6, lwd = 3)
abline(v = 0.58, col = 13, lwd = 3)
# Maximum
nlm(loglike, 0.001)$estimate      # R
#
# Example 7.22
set.seed(2)
max(runif(1000, 0, 3))
#
# Example 7.25
par(pty = "s")
par(mfrow = c(1, 2))
n <- 1000
sigma <- 1
set.seed(33)
x <- rnorm(n, 4, sigma)
mu <- seq(2, 6, length = n)
negloglikemu <- function(mu)
{ n/2 * log(2 * pi) + n/2 * log(sigma^2) + (sum(x^2)
- 2 * mu * sum(x) + n * mu^2)/(2 * sigma^2)}
EM <- nlm(negloglikemu, 2)$estimate
EM
mu1 <- 4
negloglike <- function(sigma2)
{n/2 * log(2 * pi) + n/2 * log(sigma2) +
(sum((x - mu1)^2))/(2 * sigma2)}
ES <- nlm(negloglike, 0.5)$estimate
ES
# Figure 7.6
par(mfrow=c(1,2))
plot(mu, -negloglikemu(mu), type="n")
lines(mu, -negloglikemu(mu), lwd=2)
abline(v = EM, lty = 2)
#
sigma2 <- seq(0.5, 1.5, length = 1000)
plot(sigma2, -negloglike(sigma2), type="n")
lines(sigma2, -negloglike(sigma2), lwd=2)
abline(v = ES, lty=2)
par(mfrow=c(1,1))
#
# Example 7.32
set.seed(11)
n <- 500
x <- rnorm(n, 2, 1)
mean(x)
S2u <- sum((x - mean(x))^2/n)
S2u
negloglike <- function(p)
{ (n/2)*log(2*pi) + (n/2)*log(2*p[2]) + (1/(2*p[2]))*sum((x - p[1])^2) }
nlm(negloglike, c(3, 2))$estimate
################################################################################