# Implement Fisher Scoring for linear regression

I know there is an analytic solution to the following problem (OLS). Since I try to learn and understand the principles and basics of MLE, I implemented the fisher scoring algorithm for a simple linear regression model.

$$y = X\beta + \epsilon\\ \epsilon\sim N(0,\sigma^2)$$

The loglikelihood for $\sigma^2$ and $\beta$ is given by: $$-\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln(\sigma^2)-\frac{1}{2\sigma^{2}}(y-X\beta)^{'}(y-X\beta)$$

To compute the score function $S(\theta)$, where $\theta$ is the vector of parameters $(\beta,\sigma^{2})^{'}$, I take the first partial derivatives with respect to $\beta$ and $\sigma^{2}$: \begin{align} \frac{\partial L}{\partial \beta} &= \frac{1}{\sigma^{2}}(y-X\beta)^{'}X \\[5pt] \frac{\partial L}{\partial \sigma^2} &= -\frac{N}{\sigma^{2}}+\frac{1}{2\sigma^{4}}(y-X\beta)^{'}(y-X\beta) \end{align}

Then the Fisher scoring algorithm is implemented as: $$\theta_{j+1} = \theta_{j} - (S(\theta_{j})S(\theta_{j})^{'})S(\theta_{j})$$

Please note, the following code is a very naive implementation (no stopping rule, etc.)

library(MASS)
x <- matrix(rnorm(1000), ncol = 2)
y <- 2 + x %*% c(1,3) + rnorm(500)

fisher.scoring <- function(y, x, start = runif(ncol(x)+1)){
n <- nrow(x)
p <- ncol(x)
theta <- start
score <- rep(0, p+1)
for (i in 1:1e5){
# betas
score[1:p] <- (1/theta[p+1]) * t((y - x%*%theta[1:p])) %*% x
# sigma
score[p+1] <- -(n/theta[p+1]) + (1/2*theta[p+1]^2) * crossprod(y - x %*% theta[1:p])
# new
theta <- theta - MASS::ginv(tcrossprod(score)) %*% score
}
return(theta)
}

# Gives the correct result
lm.fit(cbind(1,x), y)$coefficients # Does not give the correct result fisher.scoring(y, cbind(1,x)) # Even if you start with the correct values fisher.scoring(y, cbind(1,x), start=c(2,1,3,1))  My Question What did I miss? Where is my mistake? • A first glimpse shows$\frac{\partial L}{\partial\sigma^{2}}=-\frac{N}{{\color{red}2}\sigma^{2}}+\frac{1}{2\sigma^{4}}(y-X\beta)^{'}(y-X\beta) $. And$\theta_{j+1} = \theta_{j} {\color{red}+} (S(\theta_{j})S(\theta_{j})^{'})S(\theta_{j})$, though people would use a sum of squared first derivatives for each observation, not a squared sum. An example is here. Oct 10, 2015 at 15:44 ## 2 Answers I have fixed the code according to the suggestions by @Randel. Now it works, except$\sigma^2$takes a very long time to converge. Here is the code: fisher.scoring <- function(y, x, start = runif(ncol(x)+1)){ n <- nrow(x) p <- ncol(x) theta <- start score <- rep(0, p+1) for (i in 1:1000){ # betas score[1:p] <- (1/theta[p+1]) * t((y - x%*%theta[1:p])) %*% x # sigma score[p+1] <- -(n/(2*theta[p+1])) + (1/(2*theta[p+1]^2)) * crossprod(y - x %*% theta[1:p]) # new hessMat <- matrix(0,ncol=p+1,nrow = p+1) for(j in 1:n) { # Estimate derivative of likelihood for each observation estVec <- c((1/theta[p+1]) * t((y[j] - x[j,]%*%theta[1:p])) %*% x[j,], -(n/(2*theta[p+1])) + (1/(2*theta[p+1]^2)) * crossprod(y[j] - x[j,] %*% theta[1:p])) # Add them up as suggested to get an estimate of the Hessian. hessMat <- hessMat + estVec%*%t(estVec) } theta <- theta + MASS::ginv(hessMat) %*% score } return(theta) }  Now when I run the code I get the following: > lm.fit(cbind(1,x), y)$coefficients
x1        x2        x3
2.0136134 0.9356782 2.9666921
> fisher.scoring(y, cbind(1,x))
[,1]
[1,] 2.0136126
[2,] 0.9356782
[3,] 2.9666917
[4,] 0.6534185


I ran it again with 5000 iterations instead of 1000 and then I get:

> fisher.scoring(y, cbind(1,x))
[,1]
[1,] 2.0136133
[2,] 0.9356782
[3,] 2.9666920
[4,] 0.8962295


Hope this helps! $\sigma^2$ seems to be converging very slowly, I do not know why, I guess that is material for another question!

EDIT: Here you can see the convergence of the parameters. The number of iterations is on a log-scale. The regression coefficients take 30-40 iterations to converge, although the $\beta_1$ parameter overshoots and then comes down again, (I was not expecting to see that). The $\sigma^2$ parameter is converging rather fast as the other ones are converging, and then the convergence slows down a lot. I have no idea at the moment why this happens.

EDIT2: There was an error in the update for $\sigma^2$. See Randel's answer for the correction.

• Is there a way to solve a linear equation instead of explicitly calculating a matrix inverse in your lest line of the outer loop? I believe this is almost always a better idea when possible. Oct 11, 2015 at 15:08
• Yes @MatthewDrury you can use the solve function. I assume the code will be somewhat faster. I did not want to make too many changes in OP's code, to make the changes clear. Oct 11, 2015 at 15:15
• Fair enough. While I'm at it, thanks to both you and the poster for the extremely clear and tidy code. Oct 11, 2015 at 15:18
• @MatthewDrury you are welcome :) Oct 11, 2015 at 15:19
• My experiment shows that the estimate for $\sigma^2$ does not change much through iterations. As in your graph, if only changes 0.05, much less as compared to $\beta$. So it's almost fixed to the starting value. Oct 11, 2015 at 15:42

The short answer is that there is a bug in @guðmundur-einarsson's code.

For each observation, the score function for $\sigma^2$ is $$\frac{\partial L_i}{\partial \sigma^2} = -\frac{1}{2\sigma^{2}}+\frac{1}{2\sigma^{4}}(y_i-X_i\beta)^{'}(y_i-X_i\beta),$$ not $\frac{\partial L_i}{\partial \sigma^2} = -\frac{{\color{red}N}}{2\sigma^{2}}+\frac{1}{2\sigma^{4}}(y_i-X_i\beta)^{'}(y_i-X_i\beta)$. So just replace n in estVec with 1. Although it's not apparent, there're some clues.

• The estimate of $\sigma^2$ weirdly does not change much through iterations, although it shows some trend when ylim is automatically set to be small. This small change is because each move is divided by about $N$.
• As noted by @guðmundur-einarsson in the related question:

The convergence is rather fast in the beginning and then stagnates and continues very slowly when the regression parameters converge.

When $\beta$ does not converge (at round 4 of the x-axis in the above figures) , $\sigma^2$ depends on $\beta$, so it's changing faster as compared to that after $\beta$ converges. The log scale of iteration times seems also produce some artificial effect.

• The reason why $\beta$ can still converge with the bug can be explained by its closed-form solution does not depend on $\sigma^2$ at all.

The following figures are based on 100 iterations and the random seed is set as 1.

The lesson I learned from this is keeping debugging, although it's relatively easier in this small example.

• Thanks for pointing this out @Randel ! I deleted the other post to further reduce confusion for someone that might run into this. Oct 14, 2015 at 7:27