Projet 2

NB: Les algorithmes pourront être mis en oeuvre sous Matlab ou sous R selon votre préférance.

Catastrophes dans les mines de charbon

On voudrait modéliser l'évolution des risques d'accidents majeurs dans les mines de charbon. Pour cela, on dispose des données collectées dans le fichier coal_mine.txt. Ces données représentent la date à laquelle a eu lieu un accident impliquant plus de dix personnes dans une mine de charbon. On va supposer le modèle suivant:

Les accidents sont modélisés par un processus de Poisson.

Plus précisément, on voudrait savoir si la situation évolue en mieux au cours de l'histoire et donc on aimerait bien découvrir s'il n'y a pas plusieurs processus de Poissons successifs et des dates de changement pour passer d'un processus de Poisson au suivant. Ces dates sont appelées "breakpoints", ou "points de rupture".

Les points de ruptures sont notés t₁,...,t_d+1. On impose que le premier point de rupture est t₁=1851 et le dernier est t_d+1 = 1963, c'est à dire qu'ils sont fixés. A nous de trouver où sont les autres points de rupture ! On notera λ_i l'intensité du processus de Poisson dans l'intervalle (semi-ouvert) interval [t_i , t_i+1 ). On notera τ_j la date de chacune des 191 catastrophes. On notera n = le nombre total des catastrophes = 191 et n_i = le nombre de catastrophes dans l'intervalle [t_i , t_i+1). On prendra une loi a priori Γ (2, θ) sur les intensités des processus de Poisson, et une loi a priori Γ (2, Ψ) sur θ. On va aussi mettre une loi a priori sur les breakpoints qui soit telle que ces breakpoints se repoussent les uns les autres. On prendra la densité jointe a priori des breakpoints comme suit:

f(t2 , . . . , td ) ∝ exp(-\sum_{i=1}^d λ_i (t_i+1− t_i)) \prod_{i=1}^d λiⁿⁱ

Enfin, pour simplifier la présentation, on va collecter les paramètre dans des vecteurs:

τ = [τ₁ , · · · , τ_n ], λ = [λ₁ , · · · , λ_d ], t = [t₂ , · · · , t_d ].

Le but du projet est d'estimer θ, λ, t sachant τ. On fera un choix arbitraire mais motivé par un argument

intuitif pour Ψ.

1. Ecrire la loi jointe a posteriori de θ, λ, t sachant τ.

2. Ecrire la loi de θ sachant λ, t et τ.

3. Ecrire la loi de λ sachant θ, t et τ.

4. Construire une MCMC pour simuler θ, λ, t selon leur loi jointe a posteriori. Pour cela, on prendra une strategie de type Metropolis-Hastings (par exemple, sur θ, λ, t cycliquement pour simplifier).

5. Afficher les trajectoires de la chaine de Markov correspondant a vos simulations pour chaque parametre.

6. Afficher la fonction d'autocorrélation sur une fenêtre glissante pour chaque paramètre. Commenter.

7. Afficher un histogramme pour les lois marginales a posteriori des paramètres.

8. L'estimation est-elle très sensible au choix de Ψ ?

9. Interpréter les résultats obtenus.

10. Comment choisir le nombre de breakpoints en fin de compte ? une pénalisation de type BIC ou AIC est-elle une option défendable ?

11. Question subsidiaire: on peut aussi tenter une approche combinant Gibbs et Metropolis-hastings a partir des questions 1. et 2. Pour cela, on fera un update pour chaque paramètre successivement. θ et λ seront tirés selon leur loi a posteriori conditionnelle à t. C'est l'approche Gibbs. Les breakpoints seront perturbés selon une loi gaussienne sous contrainte de respecter l'ordre t₁ \leq ... \leq t_d+1 et le pas sera validé ou non selon l'approche Metropolis-hastings.

Ce projet est inspiré du homework2 pour le cours FMS091/MASM11 Monte Carlo-baserade statistiska metoder / Monte Carlo methods for stochastic inference de Jimmy Olsson a l'université de Lund.