# How to implement GLM computationally in C++ (or other languages)? [closed]

I want to implement the GLM model in C++ for a commercial package (ie. this is not for fun), including but not limited to normal, binomial distribution etc. I'm not so sure how the implementation should work. Say, I want to implement the linear and logistic models in a GLM framework.

I can simply implement the linear model by solving the following matrix equation for the OLS coefficients:

But I will need to iteratively estimate the MLE coefficients for a logistic model.

To me, the implementation will be something like:

if (normal in the exponential family)
{
solve the matrix equation listed above
}
else if (logistic in the exponential family)
{
iteratively solve the logistic MLE coefficients
}
else if (poisson in the ....)
{
iteratively solve the poisson MLE coefficients
}
else (.....)
{
other error distributions...
}


The GLM gives us a mathematical way of generalising the distribution models but it does not help computationally. This is like gluing distant families of different models into one model.

How does a statistical software implement a GLM model? Is there any open-source statistical library (not necessarily in C++) that shows how this can be done? Any resources (articles, papers, books) that shows how the GLM model is written?

What does the GLM model gives us how the families can be generalised computationally? Are we supposed to write an implementation for each family indepedently?

• Beware: $(X'X)^{-1}X'Y$ is the default text-book formula for computing linear regression coefficients, but it is not a good formula for computers. Nowadays people use things like QR decomposition or singular value decomposition, as these are more numerically stable. In short this is a huge project you are on, and there are many many different things that can go horribly wrong, so please reconsider whether you really should do this. Apr 14, 2015 at 7:22

While there is definitely some educational value of re-implementing GLM framework (or any other statistical framework, for that matter), I question the feasibility of this approach due to complexity and, consequently, time and efforts involved. Having said that, if you indeed want to go this route and review existing open source GLM implementations, you have, at least, the following options:

• Standard GLM implementation by R package stats. See the corresponding source code here on GitHub or by typing the function name (without parentheses) in R's command line.

• Alternative and specific GLM implementations for R include the following packages: glm2, glmnet and some others. Additionally, GLM-releated R packages are listed in this blog post.

• Excellent GLM Notes webpage (by Michael Kane and Bryan W. Lewis) offers a wealth of interesting and useful details on standard and alternative R GLM implementations aspects.

• For Julia GLM implementations, check similar to R's GLM and GLMNet packages.

• For Python GLM implementations, check the one in statsmodels library and the one in scikit-learn library (implements Ridge, OLS and Lasso - find corresponding modules).

• For .NET GLM implementations, check IMHO very interesting Accord.NET framework - the GLM source code is here on GitHub.

• For C/C++ GLM implementations, check apophenia C library (this source code seems to be relevant) and, perhaps, C++ GNU Scientific Library (GSL) (see this GitHub repo, but I was unable to find the relevant source code). Also of interest could be: this C++ IRLS GLM implementation (which uses GSL) as well as the Bayesian Object Oriented Modeling (BOOM) C++ library (GLM-focused source code is here on GitHub).