JH = function(v, w, x, y, a0, a1, a2, a3){
# x = x-coordinates
# y = y-coordinates
# a0 = alpha_0, a1 = alpha_1
E = exp(a0+a1*v+a2*w+a3*x)
L0 = -sum(E/(1+E))+sum(y)
L1 = -sum(v*E/(1+E))+sum(y*v)
L2 = -sum(w*E/(1+E))+sum(y*w)
L3 = -sum(x*E/(1+E))+sum(y*x)
L00 = -sum(E/(1+E)**2)
L01 = -sum((v*E)/(1+E)**2)
L02 = -sum((w*E)/(1+E)**2)
L03 = -sum((x*E)/(1+E)**2)
L10 = L01
L11 = -sum(v**2*E/(1+E)**2)
L12 = -sum(v*w*E/(1+E)**2)
L13 = -sum(v*x*E/(1+E)**2)
L20 = L02
L21 = L12
L22 = -sum(w**2*E/(1+E)**2)
L23 = -sum(w*x*E/(1+E)**2)
L30 = L03
L31 = L13
L32 = L23
L33 = -sum(x**2*E/(1+E)**2)
LV = c(L0, L1, L2, L3)
J = matrix(c(L00, L01, L02, L03, L10, L11, L12, L13, L20, L21, L22, L23, L30, L31, L32, L33), ncol = 4, byrow = TRUE)
list(LV = LV, Jacob = J)
}
Newton.Alg=function(v, w, x, y, ini.val, tol, maxit=100, trace = TRUE){
### initialization
err = 1
i = 1
sol.mtx = matrix(0, nrow=maxit, ncol=length(ini.val))
err.vec = rep(0, maxit)
fn.mtx = matrix(0, nrow=maxit, ncol=length(ini.val))
while(err > tol && i < maxit){
### initialization
JnFun = JH(v, w, x, y, a0 = ini.val[1],a1 = ini.val[2],a2 = ini.val[3], a3 = ini.val[4])
JMatrix = JnFun$Jacob
fn.vec = JnFun$LV
if(det(JMatrix) == 0){
cat("\n\n The Hessian matrix is singular!")
break
}
h = - solve(JMatrix)%*%fn.vec
new.val = ini.val + h
err=max(abs(h))
## store intermediate outputs
err.vec[i] = err
sol.mtx[i,] = as.vector(new.val)
fn.mtx[i,] = JH(v, w, x, y, a0 = ini.val[1],a1 = ini.val[2],a2 = ini.val[3], a3 = ini.val[4])$LV
## updating the root and the iteration ID
ini.val=new.val
i = i + 1
}
id = which(err.vec==0)[1]-1 # locate the starting rows with all zero cells
## Determinant of Hessian
SOL = sol.mtx[id,]
D = det(JH(v, w, x, y, a0 = SOL[1], a1 = SOL[2], a2 = SOL[3], a3 = SOL[4])$Jacob)
#D = det(Jacobian(sol.mtx[id,]))
fxx = (JH(v, w, x, y, a0 = SOL[1], a1 = SOL[2], a2 = SOL[3], a3 = SOL[4])$Jacob)[1,1]
##
if(trace ==TRUE){
list(solution = sol.mtx[1:id,], error = err.vec[1:id], fn.values = fn.mtx[1:id,])
} else {
if(D > 0 & fxx > 0) extreme = "local minimum"
if(D > 0 & fxx < 0) extreme = "local maximum"
if(D < 0 ) extreme = "saddle point"
list(iterations = id, solution = sol.mtx[id,],
D = D,
fxx = fxx,
extreme = extreme,
error = err.vec[id],
fn.values = fn.mtx[id,])
}
}
diabetes = read.csv("https://raw.githubusercontent.com/pengdsci/MAT325/main/w12/diabetes.csv")
diabeticStatus = diabetes$Outcome # y-variable: 1 = diabetes, 0 = no-diabetes
Glucose = diabetes$Glucose/100
BloodPressure = diabetes$BloodPressure
BMI = diabetes$BMI # x-variable (risk factor)
#glm(diabeticStatus ~ BMI, family = binomial(link = logit))
Newton.Alg(v=Glucose, w=BloodPressure, x=BMI, y=diabeticStatus, ini.val = c(-8, 4, .01, .067), tol = 10**(-8), maxit=100, trace = FALSE)
## $iterations
## [1] 5
##
## $solution
## [1] -8.80790975 4.00582201 0.01034449 0.06729461
##
## $D
## [1] 5568980656
##
## $fxx
## [1] -60.28676
##
## $extreme
## [1] "local maximum"
##
## $error
## [1] 8.817854e-09
##
## $fn.values
## [1] -1.173760e-08 6.280061e-10 -5.887523e-07 -1.458402e-07
glm(diabeticStatus ~ Glucose + BloodPressure + BMI, family = binomial(link =logit))
##
## Call: glm(formula = diabeticStatus ~ Glucose + BloodPressure + BMI,
## family = binomial(link = logit))
##
## Coefficients:
## (Intercept) Glucose BloodPressure BMI
## -8.80791 4.00582 0.01034 0.06729
##
## Degrees of Freedom: 391 Total (i.e. Null); 388 Residual
## Null Deviance: 498.1
## Residual Deviance: 371.2 AIC: 379.2
Glucose =seq(18,45, length = 50)
E = exp(-8.8 + 4*Glucose/100)
Prob.diabetes = E/(1 + E)
plot(Glucose/100, Prob.diabetes, type = "l", lwd=2, col = "blue")
BloodPressure =seq(18,45, length = 50)
E = exp(-8.8 + .01034*BloodPressure)
Prob.diabetes = E/(1 + E)
plot(BloodPressure, Prob.diabetes, type = "l", lwd=2, col = "blue")
BMI =seq(18,45, length = 50)
E = exp(-8.8 + .06729*BMI)
Prob.diabetes = E/(1 + E)
plot(BMI, Prob.diabetes, type = "l", lwd=2, col = "blue")
