# ML-NR

Jul 24th, 2020
294
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
2. # Generalized linear model: Maximum likelihood with Newton-Raphson
3.
4. ML_NR <- function(y=NULL, x=NULL, th=1e-6, start_theta=c(0,0), maxIt=100){
5.   # y: vector of  observations
6.   # x: vector of covariates
7.   # th: convergence threshold, i.e., iterative process stops
8.   #     when delta theta < th
9.   # start_theta: starting parameter
10.   # maxIt : The maximum number of iterations
11.
12.   theta_t <- start_theta
13.   it <- d.mu <- d.b <- 1
14.   theta_out <- NULL
15.   while(( abs(d.b) > th | abs(d.mu) > th ) & it < maxIt) {
16.     b <- theta_t[1]
17.     mu <- theta_t[2]
18.
19.     S <- matrix(ncol=1, nrow=2)
20.     H <- matrix(ncol=2, nrow=2)
21.     rownames(S) <- rownames(H) <- colnames(H) <- c("b", "mu")
22.
23.     Logit <- exp(b*x+mu)/(1+exp(b*x+mu))
24.     Logit2 <- exp(b*x+mu)/(1+exp(b*x+mu))^2
25.
26.     S["mu",] <- sum(y - Logit)
27.     S["b",] <- sum(x*y - Logit*x)
28.
29.     H["mu","mu"] <- -sum(Logit2)
30.     H["mu", "b"] <- -sum(Logit2*x)
31.     H["b","b"] <- -sum(Logit2*x**2)
32.     H["b", "mu"] <- -sum(Logit2*x)
33.
34.     iH <- solve(H)
35.     theta_u <- theta_t - iH%*%S
36.     d.b <- theta_u[1] - theta_t[1]
37.     d.mu <- theta_u[2] - theta_t[2]
38.
39.     theta_t <- theta_u
40.
41.     theta_out <- rbind(theta_out, c(t(theta_u),d.b, d.mu))
42.     colnames(theta_out) <- c("beta", "mu", "delta_beta", "delta_mu")
43.     it <- nrow(theta_out)
44.
45.   }
46.   theta <- theta_out
47.   vc_theta <- -iH
48.
49.   # - - - #
50.   layout(matrix(c(1,2), ncol=1, byrow=T), heights=c(5,.4))
51.   par(mar=c(5,4,4,3))
52.   plot(x=1:nrow(theta), y=theta[,3], bty="n", las=1, type="l", lwd=2, lty=5, ylab=expression(paste(Delta, "(", beta,")")), xlab="Iterations",main = "Convergence Criteria - Newton Raphson", cex.axis=.8)
53.   par(new=TRUE)
54.   plot(x=1:nrow(theta), y=theta[,4], bty="n", las=1, type="l", lwd=2, lty=3, yaxt="n", xlab="", ylab="", xaxt="n")
55.   axis(side=4, cex.axis=.8, las=1)
56.   mtext(side=4, line=2, expression(paste(Delta, "(", mu,")")))
57.
58.   par(mai=c(0,0,0,0))
59.   plot.new()
60.
61.   legend("center", legend=c(expression(paste(Delta, "(", mu,")")),
62.                             expression(paste(mu, "(", beta,")"))), lty=c(5,3), cex=.8, bty="n", horiz=TRUE)
63.   # - - - #
64.
65.   res <- list(theta=theta, "Assy variance-covariance Matrix"=vc_theta)
66.   return(res)
67. }
68.
69. file = "https://www.dropbox.com/s/f6pt59jtnrfzufi/dataAlbert2Col.txt?dl=1"
70. dat <- read.delim(file, sep=" ")
71. fit <- ML_NR(y = dat[,1], x = dat[,2])
72. fit
RAW Paste Data