SC_HW4_02
Myung
4/5/2022
Exercise J-2.4
(a)
# Data
Dose=c(0,2,10,50,250)
Test.n=c(111,105,124,104,90)
Tumor.n=c(4,4,11,13,60)
data.hw4 <- data.frame(Dose, Test.n, Tumor.n)
data <- as.matrix(data.hw4)
# E(y)=f, Jacobian, a weight matrix (inverse of variances) Function
Bino_fun <- function(beta, X){
xb = X %*% beta
f = diag(data[,2])%*%(1-exp(-xb))
J = diag(as.vector(diag(data[,2])%*%(exp(-xb))))%*%X
V = diag(as.vector(f))%*%exp(-xb)
W = diag(1/as.vector(V))
list(f=f, J=J, W=W)
}
# IRLS Function
IRLS_fun <- function(y, X, beta0, Jac, Wt = 1, maxit, tolerr, IRLS = TRUE){
fJ = function(beta, X) Jac(beta, X)
for (it in 1: maxit){
a = fJ(beta0,X)
J = a$J
f = a$f
if (IRLS==TRUE){
Wt = a$W
}
JW = t(J) %*% Wt
print(sprintf('it = %3.0f b0 = %6.6f b1 = %6.6f b2 = %6.6f, grad=%2.2e',
it,beta0[1],beta0[2],beta0[3],norm(JW %*% (y-f))))
dir = solve(JW %*% J) %*% JW %*% (y-f)
beta.n = beta0 + dir
Mre = max(abs(beta.n - beta0)/pmax(1, abs(beta.n)))
if (Mre < tolerr){ break }
beta0 = beta.n
}
return(beta.n)
}
# Start of the IRLS algorithm
X = matrix(data[,1], ncol=1)
n = dim(X)[1]
X = cbind(matrix(1,n,1),X,matrix((data[,1])^2,n,1)) # Added the constant term
# Initial Value
Ini.M = matrix(c(-1,-1,-1,-2,-10,-50,-4,-100,-2500), nrow=3, ncol=3)
R = matrix(c(log(1-4/105), log(1-11/124), log(1-13/104)), ncol=1)
Ini.V = solve(Ini.M)%*%R
Ini.V## [,1]
## [1,] 0.022934356
## [2,] 0.008191939
## [3,] -0.000119600# Run the Function
beta0 = c(0.022934356,0.008191939,-0.000119600)
y = as.matrix(data[,3])
maxit = 20
tolerr = 10^(-6)
a = IRLS_fun(y, X, beta0, Bino_fun, Wt = 1, maxit, tolerr, IRLS = TRUE)## [1] "it = 1 b0 = 0.022934 b1 = 0.008192 b2 = -0.000120, grad=5.66e+06"
## [1] "it = 2 b0 = 0.032815 b1 = 0.006930 b2 = -0.000099, grad=5.69e+06"
## [1] "it = 3 b0 = 0.034242 b1 = 0.005971 b2 = -0.000079, grad=5.78e+06"
## [1] "it = 4 b0 = 0.036136 b1 = 0.005008 b2 = -0.000059, grad=6.02e+06"
## [1] "it = 5 b0 = 0.038246 b1 = 0.004065 b2 = -0.000040, grad=6.82e+06"
## [1] "it = 6 b0 = 0.040524 b1 = 0.003173 b2 = -0.000022, grad=1.11e+07"
## [1] "it = 7 b0 = 0.042795 b1 = 0.002399 b2 = -0.000006, grad=9.73e+06"
## [1] "it = 8 b0 = 0.044651 b1 = 0.001866 b2 = 0.000005, grad=9.65e+05"
## [1] "it = 9 b0 = 0.045642 b1 = 0.001654 b2 = 0.000010, grad=8.39e+04"
## [1] "it = 10 b0 = 0.045908 b1 = 0.001628 b2 = 0.000010, grad=8.30e+02"
## [1] "it = 11 b0 = 0.045941 b1 = 0.001627 b2 = 0.000010, grad=4.89e+00"
## [1] "it = 12 b0 = 0.045944 b1 = 0.001627 b2 = 0.000010, grad=1.39e+00"# MLE for Beta0, Beta1, Beta2
a## [,1]
## [1,] 4.594453e-02
## [2,] 1.627071e-03
## [3,] 1.019261e-05(b)
# Plot
s = seq(0,250,length=250)
prob = numeric(250)
for (i in 1:250){
prob[i] = 1-exp(-a[1,]-a[2,]*s[i]-a[3,]*(s[i]^2))
}
plot(s,prob,type='l', lwd=3, col="blue", main="The Probability of Positive Response for Doses",
xlab="Doses (ppm)", ylab="Probability")
points(Dose,Tumor.n/Test.n, col="red", pch=1, cex=1)
Exercise J-2.5
(a)
# Data
data2 <- read.table("C:/RClass/blowBF.txt", header=TRUE)
# E(y)=f, Jacobian, a weight matrix (inverse of variances) Function
LR_fun <- function(beta, X){
xb = X %*% beta
f = exp(xb)/(1+exp(xb))
J = diag(as.vector(exp(xb)/(1+exp(xb))^2))%*%X
V = exp(xb)/((1+exp(xb))^2)
W = diag(1/as.vector(V))
list(f=f, J=J, W=W)
}
# IRLS Function
IRLS_fun <- function(y, X, beta0, Jac, Wt = 1, maxit, tolerr, IRLS = TRUE){
fJ = function(beta, X) Jac(beta, X)
for (it in 1: maxit){
a = fJ(beta0,X)
J = a$J
f = a$f
if (IRLS==TRUE){
Wt = a$W
}
JW = t(J) %*% Wt
print(sprintf('it = %3.0f b0 = %6.6f b1 = %6.6f b2 = %6.6f, grad=%2.2e',
it,beta0[1],beta0[2],beta0[3],norm(JW %*% (y-f))))
dir = solve(JW %*% J) %*% JW %*% (y-f)
beta.n = beta0 + dir
Mre = max(abs(beta.n - beta0)/pmax(1, abs(beta.n)))
if (Mre < tolerr){ break }
beta0 = beta.n
}
return(beta.n)
}
# Start of the IRLS algorithm
n = dim(data2)[1]
X = cbind(matrix(1,n,1), matrix(log(data2[,1]),n,1), matrix(data2[,2],n,1))
colnames(X) <- c("I","LogD", "S")
# Initial Value
model_ini <- glm(data2[,3]~log(data2[,1])+data2[,2], family="binomial", data=data2)
summary(model_ini)##
## Call:
## glm(formula = data2[, 3] ~ log(data2[, 1]) + data2[, 2], family = "binomial",
## data = data2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6629 -0.6592 -0.3509 0.5475 2.6931
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.5621 0.7499 -12.75 <2e-16 ***
## log(data2[, 1]) 3.1976 0.2999 10.66 <2e-16 ***
## data2[, 2] 4.5086 0.5159 8.74 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 856.21 on 658 degrees of freedom
## Residual deviance: 563.90 on 656 degrees of freedom
## AIC: 569.9
##
## Number of Fisher Scoring iterations: 5# Run the Function
beta0 = c(-9,3,4)
y = as.matrix(data2[,3])
maxit = 40
tolerr = 10^(-6)
a = IRLS_fun(y, X, beta0, LR_fun, Wt = 1, maxit, tolerr, IRLS = TRUE)## [1] "it = 1 b0 = -9.000000 b1 = 3.000000 b2 = 4.000000, grad=3.54e+01"
## [1] "it = 2 b0 = -9.506516 b1 = 3.177537 b2 = 4.474594, grad=1.45e+00"
## [1] "it = 3 b0 = -9.561747 b1 = 3.197443 b2 = 4.508400, grad=6.01e-03"
## [1] "it = 4 b0 = -9.562085 b1 = 3.197563 b2 = 4.508593, grad=1.49e-07"# MLE for Beta0, Beta1, Beta2
a## [,1]
## I -9.562085
## LogD 3.197564
## S 4.508593(b)
# 3D Plot
library(rgl)
options(rgl.useNULL = TRUE)
beta.MLE = matrix(c(-9.562085,3.197564,4.508593), ncol=1)
x1 = seq(0, 1, length.out=100) #S
x2 = seq(min(log(data2$D)), max(log(data2$D)), length.out=100) #log(D)
f <- function(x1, x2){
prob = exp(beta.MLE[1]+beta.MLE[2]*x2+beta.MLE[3]*x1)/
(1+exp(beta.MLE[1]+beta.MLE[2]*x2+beta.MLE[3]*x1))
prob
}
my_f <- outer(x1, x2, FUN=f)
persp3d(x1, x2, my_f, col = "blue", shade = 0.1, alpha = 0.5, xlab="S", ylab="log(D)",
zlab="prob", main="The Probability of the tree being blown down")
rglwidget(controllers = ) (c)
# Predict Probability (wind severity=0.3, tree diameter=10)
f <- function(x1, x2){
prob = exp(beta.MLE[1]+beta.MLE[2]*x2+beta.MLE[3]*x1)/
(1+exp(beta.MLE[1]+beta.MLE[2]*x2+beta.MLE[3]*x1))
prob
}
f(x1=0.3, x2=log(10))## [1] 0.3000953