<jupyter_start><jupyter_text>Load the gapminder dataWe load the data and look at the population development of a few example countries.Let's have a look at `gapminder`, which is a python dataframe<jupyter_code>gapminder.head()
gapminder.year.min(), gapminder.year.max()<jupyter_output><empty_output><jupyter_text>For a few example countries, we filter the dataframe to contain only data from that country, and plot the population<jupyter_code>for country in ['Sweden', 'Netherlands', 'United States']:
    df_country = gapminder[gapminder.country==country]
    df_country.plot(x='year', y='pop')
    plt.title(country)<jupyter_output><empty_output><jupyter_text>Exponential modelA popular model for population growth is the exponential growth model:`P(t) = P_0 e^(k(t-t_0))`We create a function that ouputs the value for an input `dt=t-t0`, and the model parameters `p0` and `k`. We use the library numpy for our calculations, so they will automatically work on arrays of data (as well as single values)<jupyter_code>import numpy as np

def model_exponential(dt, p0, k):
    return p0 * np.exp(k*dt)<jupyter_output><empty_output><jupyter_text>Let's try out this function. Suppose if we have a starting population of 1M and a growth rate of 1.01 per year. What will be the population size after 5 years?<jupyter_code>p0 = 1e6
k = 1.01
print(model_exponential(5, p0, k))<jupyter_output>156022464.48639497<jupyter_text>We can also apply the function on an array. So we calculate it for each year in one go:<jupyter_code>dt = np.arange(0,6)
print(dt)
population = model_exponential(dt, p0, k)
print(population)<jupyter_output>[1.00000000e+06 2.74560102e+06 7.53832493e+06 2.06972326e+07
 5.68263428e+07 1.56022464e+08]<jupyter_text>These numbers are a bit difficult to read, so we will plot in a barchart them using `matplotlib.pyplot`<jupyter_code>plt.bar(x=dt, height=population)
plt.show()<jupyter_output><empty_output><jupyter_text>As the numbers grow exponentially, it may be useful to have a log scale axis in our plot:<jupyter_code>plt.bar(x=dt, height=population)
plt.yscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>Now we will try to find parameters to fit the exponential curve on the growth data of the US. Let's first get the relevant data from the dataframe and convert the data to numpy arrays:<jupyter_code>df_us = gapminder[gapminder.country=='United States']
df_us.sort_values('year')
df_us
t = df_us['year'].values
dt = t - t.min()
y_observed = df_us['pop'].values

plt.scatter(t, y_observed)
plt.show()<jupyter_output><empty_output><jupyter_text>We use the first and last points to guess p0 and k. We use these guessed parameters to predict the population for each t.<jupyter_code># Initial guess for paramters:
p0_guessed = y_observed.min()
k_guessed = np.log(y_observed[-1]/p0_guessed)/dt[-1]

print('Guess for p0: {}'.format(p0_guessed))
print('Guess for k: {}'.format(k_guessed))

y_guessed = model_exponential(dt, p0_guessed, k_guessed)

plt.plot(t, y_guessed)
plt.show()<jupyter_output>Guess for p0: 157553000
Guess for k: 0.011778421587289244<jupyter_text>We can plot both the model prediction and the observed values in one plot:<jupyter_code>plt.plot(t, y_guessed, label='model', color='blue')
plt.scatter(t, y_observed, label='observed', color='red')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>ExcerciseMake an exponential model with guessed parameters for your country of choice.<jupyter_code># Fill in<jupyter_output><empty_output><jupyter_text>Now we will use the curve fitting method of scipy. We can give it our inital guesses of the parameters.<jupyter_code>from scipy.optimize import curve_fit
curve_fit?
(p0_fitted, k_fitted), pcov = curve_fit(model_exponential, dt, y_observed, p0=[p0_guessed, k_guessed])
print(p0_fitted, k_fitted)
y_fitted = model_exponential(dt, p0_fitted, k_fitted)
plt.plot(t, y_fitted, label='model', color='blue')
plt.scatter(t, y_observed, label='observed', color='red')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Logistic growth modelThe logistic growth model models an S curve, as opposed to the J curve (ever growing more!) of the exponential model. It's modeled by:`P(t) = P_m / (1 + (P_m / P_0 - 1) e^(-k(t-t0)))`Note that next to the initial population size P0, we also have the parameter P_m denoting the maximum population ExcerciseDefine a function, similar to `model_exponential`, that returns the value of the model for a given dt, P_0, P_m and k.<jupyter_code># Solution:
def model_logistic(dt, p0, pm, k):
    return pm / (1 + (pm/p0 - 1) * np.exp(-k*dt))<jupyter_output><empty_output><jupyter_text>We can make initial guesses for the parameter, but especially for k it is a bit more tricky. So let's immediately go to the curve fitting function<jupyter_code>(p0_fitted, pm_fitted, k_fitted), pcov = curve_fit(model_logistic, dt, y_observed)<jupyter_output>C:\Users\DafnevanKuppevelt\Anaconda3\lib\site-packages\scipy\optimize\minpack.py:787: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)<jupyter_text>Hm, it gives us a warning that the fitting could not be optimized. Let's help the curve fitting with some rough estimates. P_0 and P_1 are easy to guess, they are the minimum and maximum value. For k, we just start from 1.__NB:__ It is possible to work out the math, and guess the value of k from a third point.<jupyter_code>p0_guessed = y_observed.min()
pm_guessed = y_observed.max()
k_guessed = 1

(p0_fitted, pm_fitted, k_fitted), pcov = curve_fit(model_logistic, dt, y_observed, p0 = (p0_guessed, pm_guessed, k_guessed)) 
print(p0_fitted, pm_fitted, k_fitted)<jupyter_output>161519617.0014204 614264750.419242 0.01792013516982666<jupyter_text>Exercise:plot the observed population of your favorite country against those predicted by the logistic model<jupyter_code># Solution
y_fitted = model_logistic(dt, p0_fitted, pm_fitted, k_fitted)
plt.plot(t, y_fitted, label='model', color='blue')
plt.scatter(t, y_observed, label='observed', color='red')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Lognormal distribution fittingThe module `scipy.stats` offers fitting functions for different probability distributions. In this case, we don't have to define a function for the curve ourselves.We try this out on all population numbers for one year.<jupyter_code>df_2007 = gapminder[gapminder.year==2007]<jupyter_output><empty_output><jupyter_text>Let's plot the distribution of the population. It is very skewed, so we also plot the log of the population.<jupyter_code>population = df_2007['pop'].values
log_population = np.log10(population)

plt.hist(population)
plt.title('Histogram of population')
plt.show()

plt.hist(log_population)
plt.title('Histogram of log of population')
plt.show()<jupyter_output><empty_output><jupyter_text>We can also make histogram with log-scaled bins as follows:<jupyter_code>x = np.logspace(log_population.min(), log_population.max(), 30)
plt.hist(population, bins=x)
plt.xscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>Maybe we suspect that the distribution has a log-normal shape. We fit a the distribution with scipy.<jupyter_code>from scipy import stats
shape, loc, scale = stats.lognorm.fit(population)
print(shape, loc, scale)<jupyter_output>1.5395193510672078 47670.24312468737 11559300.419385288<jupyter_text>Now we can plot the probability distribution function of our fitted distribution, together with the histogram. We need to make sure that the histogram is also plotted as density.<jupyter_code>rv = stats.lognorm(shape, loc, scale)
plt.hist(population, bins=x, density=True)
plt.plot(x, rv.pdf(x))
plt.xscale('log')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10212270256277504, 'fraction_numerical': 0.03714729484856329, 'mean_word_length': 6.665674603174603, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 79, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '101874', 'license': 'Apache-2.0', 'n_tokens_mistral': 2892, 'n_tokens_neox': 2624, 'n_words': 797, 'path': 'PBL-curve-fitting.ipynb', 'repo_name': 'escience-academy/fitting-data-material', 'size': '7726'}