Parallel Model Fitting

In this notebook, we’ll fit a simple ageing model to all of the counties in the United States. As before, we’ll use scipy.optimize to perform the fitting, but we’ll use python’s multiprocessing library to perform these optimizations in parallel.

When to use this technique

This technique is appropriate when we are modeling a large number of entirely independent but structurally identical systems. In this example, we’re conceptualizing the population of counties to be influenced by aging and exogenous migration patterns. If we were to attempt to link the models together, for instance by specifying that the outmigration from one county needed to be accounted for as the inmigration to another county, we would need to perform a single large-scale optimization, or some form of hybrid.

%pylab inline
import pandas as pd
import pysd
import scipy.optimize
import multiprocessing
import numpy as np
import seaborn
Populating the interactive namespace from numpy and matplotlib



The first ingredient theat we’ll use is census data from the 2000 and 2010 census:

data = pd.read_csv('../../data/Census/Males by decade and county.csv', header=[0,1], skiprows=[2])
Unnamed: 0_level_0 Unnamed: 1_level_0 2000 2010
Unnamed: 0_level_1 Unnamed: 1_level_1 dec_1 dec_2 dec_3 dec_4 dec_5 dec_6 dec_7 dec_8 dec_9 dec_1 dec_2 dec_3 dec_4 dec_5 dec_6 dec_7 dec_8 dec_9
0 1 1 3375 3630 2461 3407 3283 2319 1637 825 284 3867 4384 3082 3598 4148 3390 2293 1353 454
1 1 3 9323 10094 7600 9725 10379 8519 6675 4711 1822 11446 12006 9976 11042 12517 12368 10623 6307 2911
2 1 5 2002 2198 2412 2465 2178 1699 1026 689 301 1673 1739 2260 2208 2233 1910 1490 739 324
3 1 7 1546 1460 1680 1762 1624 1237 774 475 187 1471 1577 1798 2016 1928 1581 1140 579 211
4 1 9 3741 3615 3393 3901 3773 3007 2227 1269 550 3741 4252 3312 3719 4129 3782 3052 1723 652


The model will be a simple ageing chain that groups individuals into 10 year cohorts.

model = pysd.read_vensim('../../models/Aging_Chain/Aging_Chain.mdl')

The Recipe

As in our other optimization problems, we’ll construct an error function that calculates the sum of squared errors between our model prediction and the measured data. We also construct a helper function called fit which basically makes the call to the optimizer and formats the result in something that we can aggregate into a Pandas DataFrame.

param_names = ['dec_%i_loss_rate'%i for i in range(1,10)]

def error(param_vals, measurements):
    predictions =, param_vals)),

    errors = predictions - measurements['2010']
    return sum(errors.values[1:]**2) #ignore first decade: no birth info

def fit(row):
    res = scipy.optimize.minimize(error, args=row,
    return pd.Series(index=['dec_%i_loss_rate'%i for i in range(1,10)], data=res['x'])

At this point, fitting the model is a simple matter of applying the fit function to the data:

county_params = data.apply(fit, axis=1)

On my 2014 era machine, this optimization takes about half an hour.

We can plot the distributions of the fit parameters for each of the counties in a histogram, to get a sense of the result. (Here we’re ignoring the first decade, which will not have reasonable parameters, as we have no information about births to the system.)

df2 = county_params.drop('dec_1_loss_rate',1)
df2.plot(kind='hist', bins=np.arange(-.15,.4,.01), alpha=.4, histtype='stepfilled')
plt.title('Fit yearly loss rates from each US county\n by age bracket from 2000 to 2010', fontsize=16)
plt.ylabel('Number of Counties', fontsize=16)
plt.xlabel('Yearly Loss Rate in 1% Brackets', fontsize=16)
plt.legend(frameon=False, fontsize=14)

Executing the optimization in parallel

We can take advantage of the multicore nature of most modern machines by using python’s multiprocessing module to distribute the various counties between each of the cores we have available for the calculation. The basic structure for this piece of code comes from this gist. We are essentially creating a helper function that will apply the fit function to a subset of the census DataFrame, and calling this function once on each of our worker nodes.


def _apply_df(args):
    df, func, kwargs = args
    return df.apply(func, **kwargs)

def apply_by_multiprocessing(df, func, workers, **kwargs):
    pool = multiprocessing.Pool(processes=workers)
    result =, [(d, func, kwargs) for d in np.array_split(df, workers)])
    return pd.concat(list(result))

county_params = apply_by_multiprocessing(data[:10], fit, axis=1, workers=4)