## Chapter 3
## November 16, 2007
## Alan T. Arnholt
#
# Example 3.8
choose(9,2)*choose(7,3)*choose(4,4)
prod(9:1)/(prod(2:1)*prod(3:1)*prod(4:1))
choose(9,2)*(prod(7:1)/(prod(3:1)*prod(4:1)))
#
# Example 3.9
for (m in seq(10,50,5))
print(c(m, 1 - prod(365:(365-m+1))/365^m))
# or
m <- seq(10,50,5)
P.E <- function(m){c(m,1-prod(365:(365-m+1))/365^m)}
t(sapply(m,P.E))
#
# Example 3.15
actual <- sample(1:3, 10000, replace = TRUE)
aguess <- sample(1:3, 10000, replace = TRUE)
equals <- (actual == aguess)
PNoSwitch <- sum(equals)/10000
not.eq <- (actual != aguess)
PSwitch <- sum(not.eq)/10000
Probs <- c(PNoSwitch, PSwitch)
names(Probs) <- c("P(Win no Switch)", "P(Win Switch)")
Probs
#
# Example 3.17
x <- 0:3
fx <- c(1/8,3/8,3/8,1/8)
Fx <- c(1/8,4/8,7/8,1) # or Fx <- cumsum(fx)
par(mfrow=c(1,2),pty="s")
plot(x, fx, type="h", xlab="x", ylab="P(X=x)",
xlim=c(0,3), ylim=c(0,.4), xaxt="n", yaxt="n")
axis(1,at=c(0,1,2,3),labels=c(0,1,2,3),las=1)
axis(2,at=c(1/8,3/8),labels=c("1/8","3/8"),las=1)
title("PDF")
plot(x, Fx,type="n", xlab="x", ylab="F(x)",
xlim=c(-1,5), ylim=c(0,1),yaxt="n")
axis(2,at=c(1/8,4/8,7/8,1),labels=c("1/8","4/8","7/8","1"),las=1)
segments(-1,0,0,0)
segments(0:4,c(Fx,1),1:5,c(Fx,1))
lines(x,Fx,type="p",pch=16)
segments(-1,1,5,1,lty=2)
title("CDF")
par(mfrow=c(1,1))
#
# Example 3.18
x <- c(1,5,30)
px <- c(0.5,0.45,0.05)
EX <- sum(x*px)
WM <- weighted.mean(x,px)
c(EX,WM)
#
# Example 3.19
x <- c(1,5,30)
px <- c(0.5,0.45,0.05)
EgX <- sum((x-5)*px)
WgM <- weighted.mean((x-5),px)
c(EgX,WgM)
#
# Example 3.20
#
# part a.
fx <- function(x){1-x^2}
k <- 1/integrate(fx,-1,1)$value
k
# part c.
fc <- function(x){.75*(1-x^2)}
integrate(fc,-.5,1)$value
# part d.
par(mfrow=c(1,2), pty="s")
x <- seq(-1,1,0.01)
y <- 3/4*(1-x^2)
plot(x, y, xlim=c(-2,2), ylim=c(0,1), type="l", xlab="x",
ylab="f(x)")
segments(-2,0,-1,0)
segments(1,0,2,0)
title("PDF for X")
y <- -x^3/4 +3*x/4+1/2
plot(x, y, xlim=c(-2,2), ylim=c(0,1), type="l", xlab="x",
ylab="F(x)")
segments(-2,0,-1,0)
segments(1,1,2,1)
title("CDF for X")
par(mfrow=c(1,1))
par(pty="m")
#
# Section 3.4.5.1 Numerical Integration with S
fx <- function(x){3/4-3/4*x^2}
integrate(fx, lower=-0.5, upper=1)
#
# Example 3.22
# part c.
curve(2*cos(2*x),0,pi/4)
abline(v=pi/12,lty=2,lwd=2)
#
# Example 3.24
x <- 0:3
px <- c(1/8,3/8,3/8,1/8)
EX <- weighted.mean(x,px)
EX2 <- weighted.mean(x^2,px)
VX <- EX2 - EX^2
sigmaX <- sqrt(VX)
c(EX,EX2,VX,sigmaX)
#
# Example 3.25
x <- 1:5
px <- x/15
EX <- sum(x*px)
sigmaX <- sqrt(sum(x^2*px) - EX^2)
X.star <- (x-EX)^3/sigmaX^3
skew <- sum(X.star*px)
skew