Statistical Computing Languages
Python
Statistical Computing Languages, RSS
[Neil D. Lawrence](http://staffwww.dcs.shef.ac.uk/people/N.Lawrence), [@lawrennd](http://twitter.com/lawrennd)
University of Sheffield
21st November 2014
In [1]:
# rerun Fernando's script to load in css
%run talktools
- Background as an Undergraduate Mechanical Engineer (1991-1994)
- Learnt BBC Basic (1983-1990), Lotus 1-2-3 (1990-1991), Quick Basic & Assembler (1991-1993), MATLAB (1994)
- Worked as Field Engineer on Oil Rigs, tried learning Borland C++ (1994-1996)
- PhD Course: Learnt C and MATLAB (1997-2000).
- Started as a lecturer in Computer Science at Sheffield 2001
- Felt obliged to learn Java
- Read Bruce Eckel's book "Thinking in Java"
- All scripts in Book were read using Python.
Python has the style of a scripting language (cf perl, bash)
In [2]:
# Open my web page and look for links (modified by code on Guido's blog)
import re
import urllib
regex = re.compile(r'href="([^"]+\.html)"')
def matcher(url, max=10):
"Print the first several URL references in a given url."
data = urllib.urlopen(url).read()
hits = regex.findall(data)
return hits
hits = matcher("http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/index.html")
for hit in hits:
print hit
- Python has the advantage that it was designed in the mid-late 1990s
- Not in the mid-late 1970s or 1980s (cf
S,MATLAB) - As a result it has a coherent and functional object system.
- It also has some attributes of functional programming.
- It is easily extended and allows for operator overloading.
- It provides an object system without imposing one (cf
java).
In [3]:
import pandas as pd
import numpy as np
class my_data(pd.DataFrame):
def pronounce(self):
print "Extending existing objects is easy, this encourages code reuse and good programming practice."
a = my_data(np.zeros((3, 3)))
a.pronounce()
a.describe()
Out[3]:
- Around turn of millenium, several researchers were looking for alternatives for
MATLAB. - Python seemed appropriate because of
- interactive shell,
- friendly scripting nature,
- well supported modern programming language
- Some, like me, realised lack of plotting, array types, numerical computing, fully functional shell was prohibitive
- Others, like Fernando Perez, John Hunter, Travis Oliphant set about adding this funcitonality.
- Targeted more at scientific computing rather than statistical computing.
- More recently Wes McKinney added
pandas: data anlaysis in python.- This gives data frame capabilities like in
R. Proving very popular for data analysis.
- This gives data frame capabilities like in
- The
scipystack catpures a range of scientific computing facilities (FFT, linear algebra etc) scikit-learnencodes models from a machine learning perspective.statsmodelsencodes models from a statistics perspective.- The Jupyter notebook and the notebook viewer provide an excellent platform for sharing analysis and code.
- For example my group makes its work available via GitHub and the nbviewer.
- Jupyter is also multilingual and has plug ins for
MATLAB,R,Julia...
- As a research group we moved from an entire
MATLABcode base to Python.- See blog post about it here
- Code is all BSD licensed (cf
Randoctavewhich are GPL licensed) - Companies deliver optimised python distributions for scientific computing (Enthought's Canopy and Continuum Analytics's Anaconda)
- I use Canopy, but we have Anaconda installed (free) on the University managed desktop in Sheffield (free). Commercial firms pay a charge for the full version.
- Like most modern languages, python is perhaps more about learning the libraries than the syntax.
- Here we use
pandasto load in data.
In [4]:
import pandas as pd # first we import the pandas library
df = pd.read_csv('../R/data/regression.csv') # now we read in the data
# gender is stored as male/female. We convert it to 0/1
df['Gender'] = df['Sex'].map( {'Female': 0, 'Male': 1} ).astype(int)
df.describe()
Out[4]:
- Plotting is done using the
matplotliblibrary.
In [5]:
# import a plotting library: matplotlib
import matplotlib.pyplot as plt
# ask for plots to appear in the notebook
%matplotlib inline
fig, ax = plt.subplots(figsize=(10, 7))
ax.scatter(df.Age, df.OI)
ax.set_xlabel('Age')
ax.set_ylabel('OI')
Out[5]:
- The plot allows us to spot the negative age value. Let's remove it.
pandasallows for removal by indexing according to 'truth values'.- These can also be used to make different colour plots.
In [6]:
df = df[df['Age']>0]
fig, ax = plt.subplots(figsize=(10, 7))
for symbol, sex in zip(['bo', 'ro'], ['Male', 'Female']):
lines= ax.semilogy(df.Age[df.Sex==sex], df.OI[df.Sex==sex], symbol,linewidth=3,markersize=10)
ax.set_xlabel('Age').set_fontsize(18)
ax.set_ylabel('OI').set_fontsize(18)
ax.set_title('Plot of age against OI on a log scale').set_fontsize(20)
First example using
statsmodels.
- Typically a lot of online resources
- Good one on regression is Connor Johnson's blog post
In [7]:
import statsmodels.formula.api as sm
df['Eins'] = 1
Y = np.log(df.OI)
Y.name = 'np.log(OI)'
X = df[['Age','Gender', 'Eins']]
result = sm.OLS( Y, X ).fit()
result.summary()
#import statsmodels.graphics.regressionplots as rp
#rp.plot_fit(results, 1)
Out[7]:
statsmodels seems designed for those with experience of statistics and R.
In [8]:
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Fit regression model (using the natural log for the dependent variable)
results = smf.ols('np.log(OI) ~ Age + Sex', data=df).fit()
# Inspect the results
results.summary()
Out[8]:
- There is even provision for loading in
Rdatasets.
In [9]:
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Load data
dat = sm.datasets.get_rdataset("Guerry", "HistData").data
# Fit regression model (using the natural log of one of the regressors)
results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()
# Inspect the results
results.summary()
Out[9]:
- In machine learning the
scikit-learnwould be more familiar.
In [10]:
from sklearn.linear_model import LinearRegression
model = LinearRegression()
X, y = df[['Age', 'Gender']], np.log(df['OI']).values[:, None]
model.fit( X, y )
score = '{0:.3f}'.format( model.score( X, y ) )
print model.coef_
- When teaching my own preference is to derive the maximum likelihood solution and have students implement it.
- Solve $$ \mathbf{X}^\top \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y} $$ for $\boldsymbol{\beta}$.
In [15]:
import numpy as np
beta = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y))
beta = pd.DataFrame(beta, index=X.columns)
beta
Out[15]:
Although Darren Wilkinson pointed out to me at the meeting, that this is more correctly done with QR decomposition (see this page for a description).
$$\mathbf{X}^\top \mathbf{X} \boldsymbol{\beta} = \mathbf{X}^\top \mathbf{y}$$$$(\mathbf{Q}\mathbf{R})^\top (\mathbf{Q}\mathbf{R})\boldsymbol{\beta} = (\mathbf{Q}\mathbf{R})^\top \mathbf{y}$$$$\mathbf{R}^\top (\mathbf{Q}^\top \mathbf{Q}) \mathbf{R} \boldsymbol{\beta} = \mathbf{R}^\top \mathbf{Q}^\top \mathbf{y}$$$$\mathbf{R}^\top \mathbf{R} \boldsymbol{\beta} = \mathbf{R}^\top \mathbf{Q}^\top \mathbf{y}$$$$\mathbf{R} \boldsymbol{\beta} = \mathbf{Q}^\top \mathbf{y}$$So if we let $\mathbf{z} = \mathbf{Q}^\top \mathbf{y}$, $$\mathbf{R} \boldsymbol{\beta} =\mathbf{z}$$
In [16]:
import numpy as np
import scipy as sp
Q, R = np.linalg.qr(X)
beta = sp.linalg.solve_triangular(R, np.dot(Q.T, y))
beta = pd.DataFrame(beta, index=X.columns)
beta
Out[16]: