whiterarebit

ML-NR

Jul 24th, 2020
294
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # author andrew.marete@canada.ca
  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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×