ols_linear_multiple_regression.c

#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multifit.h>

/* Interest rate */
double r[] = { 2.75, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.25,
               2.25, 2.25, 2,2, 2, 1.75, 1.75, 1.75, 1.75,
               1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75 };
/* Unemployment rate */
double u[] = { 5.3, 5.3, 5.3, 5.3, 5.4, 5.6, 5.5, 5.5,
               5.5, 5.6, 5.7, 5.9, 6, 5.9, 5.8, 6.1,
               6.2, 6.1, 6.1, 6.1, 5.9, 6.2, 6.2, 6.1 };
/* Stock index price */
double s[] = { 1464, 1394, 1357, 1293, 1256, 1254, 1234,
               1195, 1159, 1167, 1130, 1075, 1047, 965, 943,
               958, 971, 949, 884, 866, 876, 822, 704, 719 };

int main()
{
	int n = sizeof(r)/sizeof(double);
	gsl_matrix *X = gsl_matrix_calloc(n, 3);
	gsl_vector *Y = gsl_vector_alloc(n);
	gsl_vector *beta = gsl_vector_alloc(3);

	for (int i = 0; i < n; i++) {
		gsl_vector_set(Y, i, s[i]);

		gsl_matrix_set(X, i, 0, 1);
		gsl_matrix_set(X, i, 1, r[i]);
		gsl_matrix_set(X, i, 2, u[i]);
	}

	double chisq;
	gsl_matrix *cov = gsl_matrix_alloc(3, 3);
	gsl_multifit_linear_workspace * wspc = gsl_multifit_linear_alloc(n, 3);
	gsl_multifit_linear(X, Y, beta, cov, &chisq, wspc);

	printf("Beta:");
	for (int i = 0; i < 3; i++)
		printf("  %g", gsl_vector_get(beta, i));
	printf("\n");

	gsl_matrix_free(X);
	gsl_matrix_free(cov);
	gsl_vector_free(Y);
	gsl_vector_free(beta);
	gsl_multifit_linear_free(wspc);
}

ols_linear_multiple_regression.ipynb

ols_quadratic_fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example from: https://rosettacode.org/wiki/Multiple_regression#Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 128.81280358 -143.16202286   61.96032544]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    " \n",
    "height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,\n",
    "    1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83]\n",
    "weight = [52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93,\n",
    "    61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46]\n",
    " \n",
    "X = np.mat(height**np.arange(3)[:, None])\n",
    "y = np.mat(weight)\n",
    " \n",
    "print(y * X.T * (X*X.T).I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(height, weight)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    ,\n",
       "         1.    , 1.    , 1.    , 1.    , 1.    , 1.    , 1.    ],\n",
       "        [1.47  , 1.5   , 1.52  , 1.55  , 1.57  , 1.6   , 1.63  , 1.65  ,\n",
       "         1.68  , 1.7   , 1.73  , 1.75  , 1.78  , 1.8   , 1.83  ],\n",
       "        [2.1609, 2.25  , 2.3104, 2.4025, 2.4649, 2.56  , 2.6569, 2.7225,\n",
       "         2.8224, 2.89  , 2.9929, 3.0625, 3.1684, 3.24  , 3.3489]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[52.21, 53.12, 54.48, 55.84, 57.2 , 58.57, 59.93, 61.29, 63.11,\n",
       "         64.47, 66.28, 68.1 , 69.92, 72.19, 74.46]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}