Data Analysis and Machine Learning: Getting started, our first data and Machine Learning encounters

Introduction

Before we proceed there are several practicalities with data analysis and software tools we would like to present. These tools will help us in our understanding of various machine learning algorithms.

Our emphasis here is on understanding the mathematical aspects of different algorithms, however, where possible we will emphasize the importance of using available software. We start thus with a hands-on and top-down approach machine learning. The aim is thus to start with relevant data and use these to introduce statistical data analysis concepts and machine learning algorithms before we delve into the algorithms themselves. The examples we will use start with a simple third-order polynomial with random noise added, and using the Python software package Scikit-learn we will introduce various machine learning algorithm s to make fits of the data data and predictions. We move thereafter to more interesting cases such as the simulation of financial transactions or disease models. These are examples where we can easily set up the data and then use machine learning algorithms using included in for example scikit-learn. Another model we will consider is the so-called Ising model. Here we will use this model to produce data for selected spin configurations and attempt to classify the data. Finally, our last example consists of economic data from the OECD.

Software and needed installations

We will make intensive use of python as programming language and the myriad of available libraries. Furthermore, you will find IPython/Jupyter notebooks invaluable in your work. You can run R codes in the Jupyter/IPython notebooks, with the immediate benefit of visualizing your data.

If you have Python installed (we recommend Python3) and you feel pretty familiar with installing different packages, we recommend that you install the following Python packages via pip as

1. pip install numpy scipy matplotlib ipython scikit-learn mglearn sympy pandas pillow

For OSX users we recommend also, after having installed Xcode, to install brew. Brew allows for a seamless installation of additional software via for example

1. brew install python3

1. sudo apt-get install python3 (or python for pyhton2.7)

Python installers

1. Anaconda Anaconda is an open source distribution of the Python and R programming languages for large-scale data processing, predictive analytics, and scientific computing, that aims to simplify package management and deployment. Package versions are managed by the package management system conda
2. Enthought canopy is a Python distribution for scientific and analytic computing distribution and analysis environment, available for free and under a commercial license.

Installing R, C++, cython or Julia

You will also find it convenient to utilize R. Jupyter/Ipython notebook allows you run R code interactively in your browser. The software library R is tuned to statistically analysis and allows for an easy usage of the tools we will discuss in these texts.

To install R with Jupyter notebook following the link here

Installing R, C++, cython or Julia

For the C++ aficionados, Jupyter/IPython notebook allows you also to install C++ and run codes written in this language interactively in the browser. Since we will emphasize writing many of the algorithms yourself, you can thus opt for either Python or C++ as programming languages.

To add more entropy, cython can also be used when running your notebooks. It means that Python with the Jupyter/IPython notebook setup allows you to integrate widely popular softwares and tools for scientific computing. With its versatility, including symbolic operations, Python offers a unique computational environment. Your Jupyter/IPython notebook can easily be converted into a nicely rendered PDF file or a Latex file for further processing. For example, convert to latex as

jupyter nbconvert filename.ipynb --to latex

If you use the light mark-up language doconce you can convert a standard ascii text file into various HTML formats, ipython notebooks, latex files, pdf files etc.

Introduction to Jupyter notebook and available tools

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
import pandas as pd
from IPython.display import display
eye = np.eye(4)
print(eye)
sparse_mtx = sparse.csr_matrix(eye)
print(sparse_mtx)
x = np.linspace(-10,10,100)
y = np.sin(x)
plt.plot(x,y,marker='x')
plt.show()
data = {'Name': ["John", "Anna", "Peter", "Linda"], 'Location': ["Nairobi", "Napoli", "London", "Buenos Aires"], 'Age':[51, 21, 34, 45]}
data_pandas = pd.DataFrame(data)
display(data_pandas)

Representing data, more examples

import numpy as np
import matplotlib.pyplot as plt
from scipy import sparse
import pandas as pd
from IPython.display import display
import mglearn
import sklearn
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
x, y = mglearn.datasets.make_wave(n_samples=100)
line = np.linspace(-3,3,1000,endpoint=False).reshape(-1,1)
reg = DecisionTreeRegressor(min_samples_split=3).fit(x,y)
plt.plot(line, reg.predict(line), label="decision tree")
regline = LinearRegression().fit(x,y)
plt.plot(line, regline.predict(line), label= "Linear Regression")
plt.show()

Simple regression model

# Importing various packages
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt

x = 2*np.random.rand(100,1)
y = 4+3*x+np.random.randn(100,1)

xb = np.c_[np.ones((100,1)), x]
theta = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
xnew = np.array([[0],[2]])
xbnew = np.c_[np.ones((2,1)), xnew]
ypredict = xbnew.dot(theta)

plt.plot(xnew, ypredict, "r-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Linear Regression')
plt.show()

Simple regression model, now using scikit-learn

# Importing various packages
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

x = 2*np.random.rand(100,1)
y = 4+3*x+np.random.randn(100,1)
linreg = LinearRegression()
linreg.fit(x,y)
xnew = np.array([[0],[2]])
ypredict = linreg.predict(xnew)

plt.plot(xnew, ypredict, "r-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Random numbers ')
plt.show()

Simple regression model with gradient descent

# Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt

x = 2*np.random.rand(100,1)
y = 4+3*x+np.random.randn(100,1)

xb = np.c_[np.ones((100,1)), x]
theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
print(theta_linreg)
theta = np.random.randn(2,1)

eta = 0.1
Niterations = 1000
m = 100

for iter in range(Niterations):
    gradients = 2.0/m*xb.T.dot(xb.dot(theta)-y)
    theta -= eta*gradients

print(theta)
xnew = np.array([[0],[2]])
xbnew = np.c_[np.ones((2,1)), xnew]
ypredict = xbnew.dot(theta)
ypredict2 = xbnew.dot(theta_linreg)
plt.plot(xnew, ypredict, "r-")
plt.plot(xnew, ypredict2, "b-")
plt.plot(x, y ,'ro')
plt.axis([0,2.0,0, 15.0])
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Random numbers ')
plt.show()

Simple regression model with stochastic gradient descent

# Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import SGDRegressor

x = 2*np.random.rand(100,1)
y = 4+3*x+np.random.randn(100,1)

xb = np.c_[np.ones((100,1)), x]
theta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
print(theta_linreg)
sgdreg = SGDRegressor(n_iter = 50, penalty=None, eta0=0.1)
sgdreg.fit(x,y.ravel())
print(sgdreg.intercept_, sgdreg.coef_)

Polynomial regression

Predator-Prey model from ecology

Case study from Hudson bay

Hudson bay data

Plotting the data

import numpy as np
from  matplotlib import pyplot as plt

# Load in data file
data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
# Make arrays containing x-axis and hares and lynx populations
year = data[:,0]
hares = data[:,1]
lynx = data[:,2]

plt.plot(year, hares ,'b-+', year, lynx, 'r-o')
plt.axis([1900,1920,0, 100.0])
plt.xlabel(r'Year')
plt.ylabel(r'Numbers of hares and lynx ')
plt.legend(('Hares','Lynx'), loc='upper right')
plt.title(r'Population of hares and lynx from 1900-1920 (x1000)}')
plt.savefig('Hudson_Bay_data.pdf')
plt.savefig('Hudson_Bay_data.png')
plt.show()

Hares and lynx in Hudson bay from 1900 to 1920

Why now create a computer model for the hare and lynx populations?

The traditional (top-down) approach

Basic mathematics notation

Basic dynamics of the population of hares

Basic dynamics of the population of lynx

Evolution equations

Adapt the model to the Hudson Bay case

The program

import numpy as np
import matplotlib.pyplot as plt

def solver(m, H0, L0, dt, a, b, c, d, t0):
    """Solve the difference equations for H and L over m years
    with time step dt (measured in years."""

    num_intervals = int(m/float(dt))
    t = np.linspace(t0, t0 + m, num_intervals+1)
    H = np.zeros(t.size)
    L = np.zeros(t.size)

    print('Init:', H0, L0, dt)
    H[0] = H0
    L[0] = L0

    for n in range(0, len(t)-1):
        H[n+1] = H[n] + a*dt*H[n] - b*dt*H[n]*L[n]
        L[n+1] = L[n] + d*dt*H[n]*L[n] - c*dt*L[n]
    return H, L, t

# Load in data file
data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
# Make arrays containing x-axis and hares and lynx populations
t_e = data[:,0]
H_e = data[:,1]
L_e = data[:,2]

# Simulate using the model
H, L, t = solver(m=20, H0=34.91, L0=3.857, dt=0.1,
                 a=0.4807, b=0.02482, c=0.9272, d=0.02756,
                 t0=1900)

# Visualize simulations and data
plt.plot(t_e, H_e, 'b-+', t_e, L_e, 'r-o', t, H, 'm--', t, L, 'k--')
plt.xlabel('Year')
plt.ylabel('Numbers of hares and lynx')
plt.axis([1900, 1920, 0, 140])
plt.title(r'Population of hares and lynx 1900-1920 (x1000)')
plt.legend(('H_e', 'L_e', 'H', 'L'), loc='upper left')
plt.savefig('Hudson_Bay_sim.pdf')
plt.savefig('Hudson_Bay_sim.png')
plt.show()

The plot

If we perform a least-square fitting, we can find optimal values for the parameters \( a \), \( b \), \( d \), \( c \). The optimal parameters are \( a=0.4807 \), \( b=0.02482 \), \( d=0.9272 \) and \( c=0.02756 \). These parameters result in a slightly modified initial conditions, namely \( H(0) = 34.91 \) and \( L(0)=3.857 \). With these parameters we are now ready to solve the equations and plot these data together with the experimental values.

Linear regression in Python

import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
import sklearn
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor


data = np.loadtxt('src/Hudson_Bay.csv', delimiter=',', skiprows=1)
x = data[:,0]
y = data[:,1]
line = np.linspace(1900,1920,1000,endpoint=False).reshape(-1,1)
reg = DecisionTreeRegressor(min_samples_split=3).fit(x.reshape(-1,1),y.reshape(-1,1))
plt.plot(line, reg.predict(line), label="decision tree")
regline = LinearRegression().fit(x.reshape(-1,1),y.reshape(-1,1))
plt.plot(line, regline.predict(line), label= "Linear Regression")
plt.plot(x, y, label= "Linear Regression")
plt.show()

Linear Least squares in R

HudsonBay = read.csv("src/Hudson_Bay.csv",header=T)
fix(HudsonBay)
dim(HudsonBay)
names(HudsonBay)
plot(HudsonBay$Year, HudsonBay$Hares..x1000.)
attach(HudsonBay)
plot(Year, Hares..x1000.)
plot(Year, Hares..x1000., col="red", varwidth=T, xlab="Years", ylab="Haresx 1000")
summary(HudsonBay)
summary(Hares..x1000.)
library(MASS)
library(ISLR)
scatter.smooth(x=Year, y = Hares..x1000.)
linearMod = lm(Hares..x1000. ~ Year)
print(linearMod)
summary(linearMod)
plot(linearMod)
confint(linearMod)
predict(linearMod,data.frame(Year=c(1910,1914,1920)),interval="confidence")

Non-Linear Least squares in R

set.seed(1485)
len = 24
x = runif(len)
y = x^3+rnorm(len, 0,0.06)
ds = data.frame(x = x, y = y)
str(ds)
plot( y ~ x, main ="Known cubic with noise")
s  = seq(0,1,length =100)
lines(s, s^3, lty =2, col ="green")
m = nls(y ~ I(x^power), data = ds, start = list(power=1), trace = T)
class(m)
summary(m)
power = round(summary(m)$coefficients[1], 3)
power.se = round(summary(m)$coefficients[2], 3)
plot(y ~ x, main = "Fitted power model", sub = "Blue: fit; green: known")
s = seq(0, 1, length = 100)
lines(s, s^3, lty = 2, col = "green")
lines(s, predict(m, list(x = s)), lty = 1, col = "blue")
text(0, 0.5, paste("y =x^ (", power, " +/- ", power.se, ")", sep = ""), pos = 4)

Example: ecoli lab experiment

The program

import numpy as np

t = np.linspace(0, 10, 21)  # 20 intervals in [0, 10]
dt = t[1] - t[0]
N = np.zeros(t.size)

N[0] = 1
r = 0.5

for n in range(0, N.size-1, 1):
    N[n+1] = N[n] + r*dt*N[n]
    print 'N[%d]=%.1f' % (n+1, N[n+1])

% if FORMAT != 'ipynb':

The output

N[1]=1.2
N[2]=1.6
N[3]=2.0
N[4]=2.4
N[5]=3.1
N[6]=3.8
N[7]=4.8
N[8]=6.0
N[9]=7.5
N[10]=9.3
N[11]=11.6
N[12]=14.6
N[13]=18.2
N[14]=22.7
N[15]=28.4
N[16]=35.5
N[17]=44.4
N[18]=55.5
N[19]=69.4
N[20]=86.7

% endif

Parameter estimation

A program relevant for the biological problem

import numpy as np

# Estimate r
data = np.loadtxt('ecoli.csv', delimiter=',')
t_e = data[:,0]
N_e = data[:,1]
i = 2  # Data point (i,i+1) used to estimate r
r = (N_e[i+1] - N_e[i])/(N_e[i]*(t_e[i+1] - t_e[i]))
print 'Estimated r=%.5f' % r
# Can experiment with r values and see if the model can
# match the data better

T = 1200     # cell can divide after T sec
t_max = 5*T  # 5 generations in experiment
t = np.linspace(0, t_max, 1000)
dt = t[1] - t[0]
N = np.zeros(t.size)

N[0] = 100
for n in range(0, len(t)-1, 1):
    N[n+1] = N[n] + r*dt*N[n]

import matplotlib.pyplot as plt
plt.plot(t, N, 'r-', t_e, N_e, 'bo')
plt.xlabel('time [s]');  plt.ylabel('N')
plt.legend(['model', 'experiment'], loc='upper left')
plt.show()

Simulating financial transcations

The aim here is to simulate financial transactions among financial agents using Monte Carlo methods. The final goal is to extract a distribution of income as function of the income \( m \). From Pareto's work (V. Pareto, 1897) it is known from empirical studies that the higher end of the distribution of money follows a distribution $$ w_m\propto m^{-1-\alpha}, $$ with \( \alpha\in [1,2] \). We will here follow the analysis made by Patriarca and collaborators.

Here we will study numerically the relation between the micro-dynamic relations among financial agents and the resulting macroscopic money distribution.

We assume we have \( N \) agents that exchange money in pairs \( (i,j) \). We assume also that all agents start with the same amount of money \( m_0 > 0 \). At a given 'time step', we choose randomly a pair of agents \( (i,j) \) and let a transaction take place. This means that agent \( i \)'s money \( m_i \) changes to \( m_i' \) and similarly we have \( m_j\rightarrow m_j' \). Money is conserved during a transaction, meaning that $$ \begin{equation} m_i+m_j=m_i'+m_j'. \label{eq:conserve} \end{equation} $$ The change is done via a random reassignement (a random number) \( \epsilon \), meaning that $$ \begin{equation*} m_i' = \epsilon(m_i+m_j), \end{equation*} $$ leading to $$ \begin{equation*} m_j'= (1-\epsilon)(m_i+m_j). \end{equation*} $$ The number \( \epsilon \) is extracted from a uniform distribution. In this simple model, no agents are left with a debt, that is \( m\ge 0 \). Due to the conservation law above, one can show that the system relaxes toward an equilibrium state given by a Gibbs distribution $$ \begin{equation*} w_m=\beta \exp{(-\beta m)}, \end{equation*} $$ with $$ \begin{equation*} \beta = \frac{1}{\langle m\rangle}, \end{equation*} $$ and \( \langle m\rangle=\sum_i m_i/N=m_0 \), the average money. It means that after equilibrium has been reached that the majority of agents is left with a small number of money, while the number of richest agents, those with \( m \) larger than a specific value \( m' \), exponentially decreases with \( m' \).

We assume that we have \( N=500 \) agents. In each simulation, we need a sufficiently large number of transactions, say \( 10^7 \). Our aim is find the final equilibrium distribution \( w_m \). In order to do that we would need several runs of the above simulations, at least \( 10^3-10^4 \) runs (experiments).

Simulation of Transactions

#!/usr/bin/env python
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import random

# initialize the rng with a seed
random.seed()
# Hard coding of input parameters
Agents  = 500
MCcounts = 1000
Transactions = 100000
startMoney = 1.0
Lambda = 0.0
FinancialAgents = startMoney*np.ones(Agents)
for i in range (1, MCcounts, 1):
    for j in range (1, Transactions, 1):
        agent_i = int(Agents*random.random())
        agent_j = int(Agents*random.random())
        epsilon = random.random()
        if agent_i != agent_j:
           m1 = Lambda*FinancialAgents[agent_i] + (1-Lambda)*epsilon*(FinancialAgents[agent_i] + FinancialAgents[agent_j])
           m2 = Lambda*FinancialAgents[agent_j] + (1-Lambda)*(1-epsilon)*(FinancialAgents[agent_i] + FinancialAgents[agent_j])
           FinancialAgents[agent_i] = m1
           FinancialAgents[agent_j] = m2

# the histogram of the data
n, bins, patches = plt.hist(FinancialAgents, 50, facecolor='green')

plt.xlabel('$x$')
plt.ylabel('Distribution of wealth')
plt.title(r'Money')
plt.axis([0, 10, 0, 500])
plt.grid(True)
plt.show()

We can then change our model to allow for a saving criterion, meaning that the agents save a fraction \( \lambda \) of the money they have before the transaction is made. The final distribution will then no longer be given by Gibbs distribution. It could also include a taxation on financial transactions.

The conservation law of Eq. \eqref{eq:conserve} holds, but the money to be shared in a transaction between agent \( i \) and agent \( j \) is now \( (1-\lambda)(m_i+m_j) \). This means that we have $$ \begin{equation*} m_i' = \lambda m_i+\epsilon(1-\lambda)(m_i+m_j), \end{equation*} $$ and $$ \begin{equation*} m_j' = \lambda m_j+(1-\epsilon)(1-\lambda)(m_i+m_j), \end{equation*} $$ which can be written as $$ \begin{equation*} m_i'=m_i+\delta m \end{equation*} $$ and $$ \begin{equation*} m_j'=m_j-\delta m, \end{equation*} $$ with $$ \begin{equation*} \delta m=(1-\lambda)(\epsilon m_j-(1-\epsilon)m_i), \end{equation*} $$ showing how money is conserved during a transaction. Select values of \( \lambda =0.25,0.5 \) and \( \lambda=0.9 \) and try to extract the corresponding equilibrium distributions and compare these with the Gibbs distribution. Comment your results. Extract a parametrization of the above curves, see for example Patriarca and collaborators and see if you can parametrize the high-end tails of the distributions in terms of power laws. Comment your results.

In the studies above the agents were selected randomly, irrespective of whether we allowed for saving or not during a transaction. What is often observed is that various agents tend to make preferences for for whom to interact with. We will now study the evolution of the distribution of wealth \( w_m \) by assuming that there is a likelihood $$ p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}, $$ for an interaction between agents \( i \) and \( j \) with respective wealths \( m_i \) and \( m_j \). The parameter \( \alpha > 0 \). For \( \alpha=0 \) we recover our model from part 5a). Perform the same analysis as previously with \( N=500 \) as well as with \( N=1000 \) agents and study the distribution of wealth for \( \alpha =0.5 \), \( \alpha =1.0 \), \( \alpha =1.5 \) and \( \alpha =2.0 \). You should try to reproduce Figure 1 of Goswami and Sen. Extract the tail of the distribution and see if it follows a Pareto distribution $$ w_m\propto m^{-1-\alpha}. $$ What happens if \( \alpha \gg 1 \)?

Perform the analysis with and without a saving \( \lambda \) on each transaction and comment your results. We add to the previous probability the possibility that two agents who interact have performed similar transactions earlier. That is, in addition to being financially close, we assume that the likelihood for interacting increases if two agents have interacted earlier. We add this feature by modifying the previous likelihood to $$ p_{ij} \propto \vert m_i-m_j\vert^{-\alpha}\left(c_{ij}+1\right)^{\gamma}, $$ where \( c_{ij} \) represents the number of previous interactions that have taken place between \( i \) and \( j \). The factor \( 1 \) is added in order to ensure that if they have not interacted earlier they can still interact. Perform similar studies as above with \( N=1000 \), \( \alpha=1.0 \) and \( \alpha=2.0 \) using \( \gamma = 0.0, 1.0, 2.0, 3.0 \) and \( 4.0 \). Plot the wealth distributions for these cases and try to extract eventual power law tails with and without a saving \( \lambda \) in each transaction. Comment your results and compare them with figures 5 and 6 of Goswami and Sen.

Particle in one dimension an velocity distribution

# Program to test the Metropolis algorithm with one particle at given temp in one dimension
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import random
from math import sqrt, exp, log
# initialize the rng with a seed
random.seed()
# Hard coding of input parameters
MCcycles = 100000
Temperature = 2.0
beta = 1./Temperature
InitialVelocity = -2.0
CurrentVelocity = InitialVelocity
Energy = 0.5*InitialVelocity*InitialVelocity
VelocityRange = 10*sqrt(Temperature)
VelocityStep = 2*VelocityRange/10.
AverageEnergy = Energy
AverageEnergy2 = Energy*Energy
VelocityValues = np.zeros(MCcycles)
# The Monte Carlo sampling with Metropolis starts here
for i in range (1, MCcycles, 1):
    TrialVelocity = CurrentVelocity + (2.0*random.random() - 1.0)*VelocityStep
    EnergyChange = 0.5*(TrialVelocity*TrialVelocity -CurrentVelocity*CurrentVelocity);
    if random.random() <= exp(-beta*EnergyChange):
        CurrentVelocity = TrialVelocity
        Energy += EnergyChange
        VelocityValues[i] = CurrentVelocity
    AverageEnergy += Energy
    AverageEnergy2 += Energy*Energy
#Final averages
AverageEnergy = AverageEnergy/MCcycles
AverageEnergy2 = AverageEnergy2/MCcycles
Variance = AverageEnergy2 - AverageEnergy*AverageEnergy
print(AverageEnergy, Variance)
n, bins, patches = plt.hist(VelocityValues, 400, facecolor='green')

plt.xlabel('$v$')
plt.ylabel('Velocity distribution P(v)')
plt.title(r'Velocity histogram at $k_BT=2$')
plt.axis([-5, 5, 0, 600])
plt.grid(True)
plt.show()