bwlewis / glm Goto Github PK

We'll roll all these methods in to an R package.

But, the SVD Newton method will be the primary non-distributed method. We should add the following to it:

1. R-like subset selection
2. Full set of output stats

Also note that the termination criterion is not the same as R, we should change that.

I always found your minimalist IRLS GLM algorithm at [http://bwlewis.github.io/GLM/ ] very handy to explain to students how GLMs work. I was wondering if, by any chance, you would also happen to have an equally minimalist implementation of multinomial & negative binomial GLMs available, ideally also using IRLS (as opposed to Newton-Raphson)?

These notes are incredibly useful.

I think I've found a bug though. Check "Algorithm IRLS: Iteratively reweighted least squares estimation" in the "Input" section. I think the line:

Let x_{j+1}=x_j+(A^TWA)^{−1}A^TWz.

should be

Let x_{j+1}= (A^TWA)^{−1}A^TWz.

I noticed that all algo's listed are significantly slower than the inbuilt glm.fit. Why is this? Does glm.fit use better starting values perhaps? Or is this due to the use of .Call(C_Cdqrls...) in glm.fit as opposed to the use of solve() in the irls code?

In the sparse and incremental variants, the resulting linear system is solved by Cholesky decomposition with pivoting. So $X^T X \beta = X^T y$ is solved essentially with

U <- chol(crossprod(X),pivot=TRUE)
p <- attr(U,"pivot")
backsolve(U,backsolve(U,crossprod(X,y)[p],transpose=TRUE))[p]

But shouldn't the final line be

backsolve(U,backsolve(U,crossprod(X,y)[p],transpose=TRUE))[order(p)]

to invert the original permutation?