transition starting

abc-ages/Readme.md

+# Population Model
+
+Model of the dynamics of population for French farmers that includes two
+  processes, an installation process and a retirement process.
+The population is represented by an array of age class amounts.
+The installation process is modeled by a beta distribution with 4 parameters:
+  amplitude, drop_x, mu and v. The retirement process is a linear function
+  with a threshold.
+
+See documentation in the notebook "ABC Calibration Yearly Model.ipynb"
+
+## Usage
+
+First, you need to calibrate the model with the script abc.R. You can adjust the number
+of samples `nb_simul`, then run:
+```
+Rscript abc.R
+```
+
+Then you can use the Python API exposed by `population_model.py` for using the calibrated parameters and computing the yearly dynamics.

abc-ages/abc.R

+library(EasyABC)
+library(parallel)
+library(lhs)
+library(mnormt)
+
+# This script implements the population model and calibrates its parameters by
+# using ABC Lenormand method.
+# The output will be generated in ../out
+dir.create("../out")
+
+# MSA 2008
+ages_source = c(5, 47, 194, 529, 1056, 1760, 2477, 3184, 4013, 4617, 5210, 5495, 5981, 6354, 6949, 7646, 8706, 9973, 10766, 11644, 12479, 13180, 13736, 14781, 15885, 16682, 17803, 17706, 17987, 18386, 18158, 17949, 17533, 18011, 17308, 17259, 17538, 18032, 18606, 17364, 16668, 15279, 13731, 7604, 5294, 3282, 2669, 2286, 1614, 1236, 1007, 1074, 978, 829, 789, 655, 697, 591, 496, 471, 403, 348, 323, 269, 248, 240, 220, 171, 178, 201, 248, 130, 58, 63, 43, 42, 55, 42, 30, 30, 9, 10, 15)
+# MSA 2018
+ages_target = c(1, 76, 268, 579, 1079, 1493, 1973, 2432, 3059, 3592, 4236, 4717, 5440, 5904, 6419, 6945, 7574, 7847, 8414, 8704, 8992, 8782, 8821, 9069, 9236, 9741, 10569, 11458, 12058, 12623, 12927, 13351, 13637, 14446, 15461, 16145, 17113, 17017, 17163, 17660, 17350, 17190, 16147, 12294, 10037, 6746, 5509, 4670, 3720, 2830, 2401, 2181, 1856, 1584, 1209, 790, 665, 595, 496, 390, 333, 346, 303, 246, 238, 205, 201, 175, 146, 125, 92, 82, 75, 55, 50, 36, 36, 26, 25, 18, 13, 10, 11)
+D=10 # Years delta between the two sets
+
+simulator_yearly = function(ages, x) {
+  # computation of the new age at index considering parameters "x"
+  age_process = function(index) {
+    set.seed(x[1])
+    # installation
+    amplitude = x[2]
+    drop_x=x[3]
+    alpha = x[4]*x[5]
+    beta = (1 - x[4])*x[5]
+    # retirement
+    retirement_threshold = x[6]
+    retirement_factor = x[7]
+    if (index < retirement_threshold) retirement_factor = 0
+    new_ages[index] + amplitude*dbeta(index/drop_x, alpha, beta) - new_ages[index]*retirement_factor
+  }
+  new_ages = c(0, ages[1:(length(ages)-1)])
+  # Compute all the new ages
+  unlist(lapply(seq_along(new_ages), age_process))
+  # 'unlist' is required because EasyABC requires a vector as output instead of a list
+}
+simulator = function(x) {
+    ages = ages_source
+    for (i in seq(1,D)) {
+        ages = simulator_yearly(ages, x)
+    }
+    ages
+}
+
+paramsNames = c("inst_amplitude","inst_drop","inst_beta1",
+        "inst_beta2","retire_threshold","retire_factor")
+nb_simul = 10000
+p_acc_min=0.4
+file_prefix=paste("../out/abc_",format(nb_simul),".Rdata",sep="")
+prior_unif = list(c("unif",200,400), c("unif",25,40), c("unif",0.3,.45),
+                  c("unif",4,10), c("unif",40,55), c("unif",0,0.3))
+ABC_Lenormand <- ABC_sequential(method="Lenormand", model=simulator,
+      prior=prior_unif, nb_simul=nb_simul, summary_stat_target=ages_target,
+      p_acc_min=p_acc_min,use_seed=TRUE, progress_bar=TRUE)
+# Save data
+write.csv(ABC_Lenormand$param,paste(file_prefix,"_param.csv",sep=""),row.names=F)
+write.csv(ABC_Lenormand$weights,paste(file_prefix,"_weights.csv",sep=""),row.names=F)
+# Plot results againts data
+simulated = ABC_Lenormand$stats
+png(paste(file_prefix,".png",sep=""))
+boxplot(simulated,outpch=NA)
+points(ages_source,pch=3,col="green")
+points(ages_target,pch=3,col="red")
+legend("topleft", legend=c(
+      paste("ABC Lenormand \n",format(nb_simul),
+      " simulations pacc=",format(p_acc_min),sep=""
+    ),"Source","Target"),
+    col=c("black","green","red"),cex=.8,pch=3)
+dev.off()
+# Plot posteriors
+png(paste(file_prefix,"_posteriors.png",sep=""))
+par(mfrow=c(3,2))
+for (c in 1:ncol(ABC_Lenormand$param)) {
+  plot(density(ABC_Lenormand$param[,c],weights = ABC_Lenormand$weights),xlab=paramsNames[c],main=paramsNames[c])
+}
+dev.off()

abc-ages/population_model.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+import pandas as pd
+import numpy as np
+import scipy
+from scipy import stats
+
+class PopulationModel():
+    '''
+    Model of the dynamics of population for French farmers that includes two
+     processes, an installation process and a retirement process.
+    The population is represented by an array of age class amounts.
+    The installation process is modeled by a beta distribution with 4 parameters:
+     amplitude, drop_x, mu and v. The retirement process is a linear function
+     with a threshold.
+    See documentation in the notebook "ABC Calibration Daily Model.ipynb".
+    '''
+    def __init__(self, ageclasses, params_filename, weights_filename, nbsamples=1):
+        '''
+        Initialize an instance of the model
+        :param ageclasses inital data of age classes amounts
+        :param params_filename filename of the posteriors of the calibrated
+        parameters
+        :param weights_filename filename of the weights for this posteriors
+        :param nbsamples number of samples that will be handled by this instance
+        '''
+        self.ageclasses = np.repeat([ageclasses], nbsamples, axis=0)
+        params = pd.read_csv(params_filename)
+        weights = pd.read_csv(weights_filename).iloc[:,0]
+        kernel = stats.gaussian_kde(params.transpose(), weights=weights)
+        self.parameters = kernel.resample(nbsamples).transpose()
+
+    def step_year(self):
+        '''
+        Compute the evolution of the age classes amounts for one year. The
+        new amounts will be stored in the instance, ready for a next step, and
+        will also be returned.
+        '''
+        # First, we shift the amounts by one year, and affect 0 for the
+        # new ageclass
+        self.ageclasses = np.roll(self.ageclasses,1)
+        self.ageclasses[:,0] = 0
+        # Now, we process every parameters set sampled
+        for i in range(len(self.ageclasses)):
+            [amplitude,drop_x,mu,v,r_threshold,r_factor] = self.parameters[i,:]
+            alpha = mu * v
+            beta = (1 - mu) * v
+            newages = self.ageclasses[i,:]
+            newages = map(lambda elt:
+                      elt[1]
+                      + amplitude*stats.beta.pdf((1+elt[0])/drop_x, alpha, beta)
+                       - 0 if ((1+elt[0]) < r_threshold) else elt[1]*r_factor
+                      , enumerate(newages))
+            self.ageclasses[i,:] = list(newages)
+        return self.ageclasses
+
+if __name__ == '__main__':
+    ages_msa_2008 = [5, 47, 194, 529, 1056, 1760, 2477, 3184, 4013, 4617, 5210, 5495, 5981, 6354, 6949, 7646, 8706, 9973, 10766, 11644, 12479, 13180, 13736, 14781, 15885, 16682, 17803, 17706, 17987, 18386, 18158, 17949, 17533, 18011, 17308, 17259, 17538, 18032, 18606, 17364, 16668, 15279, 13731, 7604, 5294, 3282, 2669, 2286, 1614, 1236, 1007, 1074, 978, 829, 789, 655, 697, 591, 496, 471, 403, 348, 323, 269, 248, 240, 220, 171, 178, 201, 248, 130, 58, 63, 43, 42, 55, 42, 30, 30, 9, 10, 15]
+    model = PopulationModel(ages_msa_2008,"../out/abc_10000.Rdata_param.csv", "../out/abc_10000.Rdata_weights.csv",nbsamples=1)
+    for i in range(10):
+        target = model.step_year()
+    data=pd.DataFrame(target)
+    # print(data.head())
+    data.to_csv("../out/data.csv")

ardoise.py

+import pandas as pd
+from numpy import random
+# Build cars DataFrame
+names = ['United States', 'Australia', 'Japan', 'India', 'Russia', 'Morocco', 'Egypt']
+dr =  [True, False, False, False, True, True, True]
+cpc = [809, 731, 588, 18, 200, 70, 45]
+lst = [random.randint(10, size=random.randint(1,5)) for i in range(len(names))]
+dict = { 'country':names, 'cars_per_cap':cpc, "lst": lst}
+cars = pd.DataFrame(dict)
+# print(cars)
+
+# Definition of row_labels
+row_labels = ['US', 'AUS', 'JAP', 'IN', 'RU', 'MOR', 'EG']
+
+# Specify row labels of cars
+cars.index = row_labels
+
+# Print cars again
+cars_sample = cars.sample(3)
+cars_sample["country"] = "foo"
+cars_sample["lst"] = cars_sample["lst"].apply(lambda x: [5 for i in range(len(x))])
+print(cars)
+print("---")
+print(cars_sample)
+cars = cars.drop(list(cars_sample.index.values))
+print(cars)
+print("---")
+cars = pd.concat([cars_sample, cars], axis=0)
+print(cars)
\ No newline at end of file

src/transition.py

+import pandas as pd
+import numpy as np
+
+# Here will go the transition function described in transition_age.ipynb
\ No newline at end of file