More on Carbon Taxes

I posted earlier about why carbon taxes are the way to go.

NPR just released an article that says pretty much the same thing. Carbon tax = reduced emissions and climate policy + increased revenue for the government. It also allows market-based solutions to climate change, rather than a regulatory approach.

Evaluating Optimization Algorithms in MATLAB, Python, and R

As those of you who read my last post know, I’m at the NIMBioS-CAMBAM workshop on linking mathematical models to biological data here at UT Knoxville. Day 1 (today) was on parameter estimation and model identifiability. Specifically, we (quickly) covered maximum likelihood estimation and how to program our own estimation procedures in MATLAB.

If you’ve read this blog in the past, you know I’m an open-source kind of person. Naturally, I re-programmed all of the MATLAB examples in R. When I did, I noticed something funny. The optimization procedures in MATLAB gave different estimates than those in R. I asked a post-doc there, who seemed equally stumped but did mention that R’s optimization procedures are little funky. So, I took the next logical step and programmed the ML optimization routine into Python, using Scipy and Numpy, just to double check.

The model is an SIR epidemiological model that predicts the number of Susceptible, Infected, and Recovering people with, in this case, cholera. It relies on four parameters, Bi, Bw, e, and k. I won’t give the model here, you’ll see the formula in the code below. When optimizing, I made sure that MATLAB, R, and Python all used Nelder-Mead algorithms and, when possible, equivalent ODE solvers (ode45 in MATLAB and R).

I won’t post the MATLAB code here, because I didn’t write it and it’s multiple files etc etc, but I’ve gone over it line by line to make sure it’s identical to my R and Python code. It is.
MATLAB outputs estimates for Bi (0.2896), Bw (1.0629), e (0.0066) and k (0.0001). If you want the MATLAB files, or the data, I’ll send them to you.

In R, first import the data (see the Python code below for the actual data). Then, run the ode solver and optimization code:

## Main code for the SIR model example
library(deSolve) # Load the library to solve the ode

## Set initial parameter values
Bi <- 0.75
Bw <- 0.75
e <- 0.01
k <- 1/89193.18

## Combine parameters into a vector
params <- c(Bi, Bw, e, k)
names(params) <- c('Bi', 'Bw', 'e', 'k')

# Make a function for running the ODE
SIRode <- function(t, x, params){
S <- x[1]
I <- x[2]
W <- x[3]
R <- x[4]

Bi <- params[1]
Bw <- params[2]
e <- params[3]
k <- params[4]

dS <- -Bi*S*I - Bw*S*W
dI <- Bi*S*I + Bw*S*W - 0.25*I
dW <- e*(I - W)
dR <- 0.25*I

output <- c(dS, dI, dW, dR)
list(output)
}

# Set initial conditions
I0 <- data[1]*k
R0 <- 0
S0 <- (1-I0)
W0 <- 0

initCond <- c(S0, I0, W0, R0)

# Simulate the model using our initial parameter guesses
initSim <- ode(initCond, tspan, SIRode, params, method='ode45')

plot(tspan, initSim[,3]/k, type='l')
points(tspan, data)

# Make a function for optimzation of the parameters
LLode <- function(params){

k <- params[4]

I0 <- data[1]*k
R0 <- 0
S0 <- 1 - I0
W0 <- 0

initCond <- c(S0, I0, W0, R0)

# Run the ODE
odeOut <- ode(initCond, tspan, SIRode, params, method='ode45')

# Measurement variable
y <- odeOut[,3]/k

diff <- (y - data)
LL <- t(diff) %*% diff

return(LL)
}

# Run the optimization procedure
MLresults <- optim(params, LLode, method='Nelder-Mead')

## Resimulate the ODE
estParms <- MLresults$par
estParms
MLresults$value

I0est <- data[1]*estParms[4]
S0est <- 1 - I0est
R0 <- 0
W0 <- 0
initCond <- c(S0est, I0est, W0, R0)

odeEstOut <- ode(initCond, tspan, SIRode, estParms, method='ode45')
estY <- odeEstOut[,3]/estParms[4]

plot(tspan, data, pch=16, xlab='Time', ylab='Number Infected')
lines(tspan, estY)

Running this code gives me estimates of Bi (0.30), Bw (1.05), e (0.0057), and k (0.0001).

R_OptimizeI rounded the values, but the unrounded values are even closer to MATLAB. So far, I’m fairly satisfied.

We can run the same procedure in Python:


from scipy.optimize import minimize
from scipy import integrate
import pylab as py

## load data

tspan = np.array([0, 7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84, 91, 98, 105, 112, 119, 126, 133, 140, 147, 154, 161])
data = np.array([ 113, 60, 70, 140, 385, 2900, 4600, 5400, 5300, 6350, 5350, 4400, 3570, 2300, 1900, 2200, 1700, 1170, 830, 750, 770, 520, 550, 380 ])

## define ODE equations
def SIRode(N, t, Bi, Bw, e, k):
    return(
    -Bi*N[0]*N[1] - Bw*N[0]*N[2] ,
    +Bi*N[0]*N[1] + Bw*N[0]*N[2] - 0.25*N[1] ,
    +e*(N[1] - N[2]) ,
    +0.25*N[1]
    )

# Set parameter values
Bi = 0.75
Bw = 0.75
e = 0.01
k = 1/89193.18

# Set initial conditions
I0 = data[0]*k
S0 = 1 - I0
N = [S0, I0, 0]

# Run the ode
Nt = integrate.odeint(SIRode, N, tspan, args=(Bi, Bw, e, k))

# Get the second column of data corresponding to I
INt = [row[1] for row in Nt]
INt = np.divide(INt, k)

py.clf()
py.plot(tspan, data, 'o')
py.plot(tspan, INt)
py.show()

def LLode(x):
    Bi = x[0]
    Bw = x[1]
    e = x[2]
    k = x[3]

    I0 = data[0]*k
    S0 = 1-I0
    N0 = [S0, I0, k]

    Nt = integrate.odeint(SIRode, N0, tspan, args=(Bi, Bw, e, k))

    INt = [row[1] for row in Nt]
    INt = np.divide(INt, k)

    difference = data - INt

    LL = np.dot(difference, difference)

    return LL

x0 = [Bi, Bw, e, k]

results = minimize(LLode, x0, method='nelder-mead')

print results.x

estParams = results.x

Bi = estParams[0]
Bw = estParams[1]
e = estParams[2]
k = estParams[3]

I0 = data[0]*k
S0 = 1 - I0
N = [S0, I0, 0]

Nt = integrate.odeint(SIRode, N, tspan, args=(Bi, Bw, e, k))

INt = [row[1] for row in Nt]
INt = np.divide(INt, k)

py.clf()
py.plot(tspan, data, 'o')
py.plot(tspan, INt)
py.show()

Python gave me estimates of Bi (0.297), Bw (1.106), e (0.0057), and k (0.00001).

Python_optimizeOverall, I have to say I’m pretty satisfied with the performance of both R and Python. I also didn’t find programming these sorts of optimization procedures in R or Python to be any more difficult than in MATLAB (discounting for the fact that I’m not terribly familiar with MATLAB, but I’m also only somewhat familiar with Python, so they’re roughly equivalent here). I initially wrote this post because I could not get R to sync up with MATLAB, so I tried it in Python and got the same results as in R. I then found out that something was wrong with the MATLAB code, so that all three matched up pretty well.

I’ll say it again, R + Python = awesome.

UPDATE:

I’ve had a couple of requests for time comparisons between MATLAB, R, and Python. I don’t have those for MATLAB, I don’t know MATLAB well.

For Python, if I wrap the mimize function in:

t0 = time.time()
miminize....# run the optimizer
t1 = time.time()
print t1 - t0

I get 3.17 seconds.

In R, if I use system.time( ) to time the optim( ) function, I get about 39 seconds. That pretty much matches my feeling that R is just laboriously slow compared with how quickly Python evaluates the function.

Happy Monday From NIMBioS

Happy Monday from UT Knoxville! I start the NIMBioS – CAMBAM workshop on mathematical models in biology today and it runs for two weeks. I’m excited and concerned. I’m excited because I stand to learn a lot. I’m concerned because I just got my experiments up and running on Saturday before I left.

130615_001

There is a distinct lack of heaters and power….

The seeds are planted, hopefully growing, but there’s still no power in the garden to run my heaters. I’m stressing because I’m out of town on these crucial two weeks: Will I get power or will I need to redo¬†everything? Will my seeds germinate? Did I screw i all up? I’m out of town and can’t watch it or fix it, so… bonne chance.

130615_004

Grow, baby.. grow

So, to make myself feel better about it, here’s a picture of the Happy Monday Dog.

IMG958860

She is also stressed… Can you tell?

Happy Monday Dog!

I did have a long post written about the NSA phone scandal and why it doesn’t bother me, but I had spent 10 hours weighing plant seeds in the lab yesterday (and I’m only halfway done) and it was a little crazy. My friend convinced me to not post it, so instead I leave you with the Happy Monday Dog, because she makes everything seem better.

Government Watchdog

A Government Watchdog

More seed weighing today!