Fisher Scoring Algorithm (R version)

Fisher_Scoring.R

getCoefficients <- function(design_matrix, response_vector, epsilon=.0001) {
    # =========================================================================
    # design_matrix      `X`     => n-by-(p+1)                                |
    # response_vector    `y`     => n-by-1                                    |
    # probability_vector `p`     => n-by-1                                    |
    # weights_matrix     `W`     => n-by-n                                    |
    # epsilon                    => threshold above which iteration continues |
    # =========================================================================
    # n                          => # of observations                         |
    # (p + 1)                    => # of parameterss, +1 for intercept term   |
    # =========================================================================
    # U => First derivative of Log-Likelihood with respect to                 |
    #      each beta_i, i.e. `Score Function`: X_transpose * (y - p)          |
    #                                                                         |
    # I => Second derivative of Log-Likelihood with respect to                |
    #      each beta_i. The `Information Matrix`: (X_transpose * W * X)       |
    #                                                                         |
    # X^T*W*X results in a (p+1)-by-(p+1) matrix                              |
    # X^T(y - p) results in a (p+1)-by-1 matrix                               |
    # (X^T*W*X)^-1 * X^T(y - p) results in a (p+1)-by-1 matrix                |
    # ========================================================================|
    X <- as.matrix(design_matrix)
    y <- as.matrix(response_vector)
    
    # initialize logistic function used for Scoring calculations =>
    pi_i <- function(v) return(exp(v)/(1 + exp(v)))
    
    # initialize beta_0, p_0, W_0, I_0 & U_0 =>
    beta_0 <- matrix(rep(0, ncol(X)), nrow=ncol(X), ncol=1, byrow=FALSE, dimnames=NULL)
    p_0    <- pi_i(X %*% beta_0)
    W_0    <- diag(as.vector(p_0*(1-p_0)))
    I_0    <- t(X) %*% W_0 %*% X
    U_0    <- t(X) %*% (y - p_0)
    
    # initialize variables for iteration =>
    beta_old                   <- beta_0
    iter_I                     <- I_0
    iter_U                     <- U_0
    iter_p                     <- p_0
    iter_W                     <- W_0
    fisher_scoring_iterations  <- 0
    
    # iterate until difference between abs(beta_new - beta_old) < epsilon =>
    while(TRUE) {
        
        # Fisher Scoring Update Step =>
        fisher_scoring_iterations <- fisher_scoring_iterations + 1
        beta_new <- beta_old + solve(iter_I) %*% iter_U
        
        if (all(abs(beta_new - beta_old) < epsilon)) {
            model_parameters  <- beta_new
            fitted_values     <- pi_i(X %*% model_parameters)
            covariance_matrix <- solve(iter_I)
            break
            
        } else {
            iter_p   <- pi_i(X %*% beta_new)
            iter_W   <- diag(as.vector(iter_p*(1-iter_p)))
            iter_I   <- t(X) %*% iter_W %*% X
            iter_U   <- t(X) %*% (y - iter_p)
            beta_old <- beta_new
        }
    }
    
    summaryList <- list(
        'model_parameters'=model_parameters,
        'covariance_matrix'=covariance_matrix,
        'fitted_values'=fitted_values,
        'number_iterations'=fisher_scoring_iterations
    )
    return(summaryList)
}