TP2 : Estimation des paramètres de Malkmus à partir de statistiques de k

Ce TP fait suite au TP1 où le modèle de Malkmus est ajusté aux données de transmissivité générées par le programme sln-slab. On va ici essayer d'estimer les paramètres du modèle de Malkmus directement à partir de statistiques sur les coefficients d'absorption k, à l'aide du programme sln-stat.

On suppose que l'on travaille dans un répertoire ${HOME}/mon_nom/tp_star-line/TP2 conformément aux instructions données dans la présentation générale du TP.

export LANG=en_US.UTF-8

Partie 1 : Générer des statistiques de k

Le programme sln-stat est livré avec la bibliothèque star-line. Il implémente un calcul par Monte Carlo de k_barre et k^2_barre, définies comme les moyennes de k_nu et de k_nu^2 sur une bande étroite. Ces moyennes peuvent être écrites comme les espérances de variables aléatoires E[K_nu] et E[K_nu^2] (les lettres capitales étant réservées aux variables aléatoires). On a donc que les K_nu sont dont distribuées. Les distributions de K, ou k-distributions, sont un concept au coeur de ce TP. Une façon standard d'étudier les distributions de la variable aléatoire K_nu est de commencer par calculer un spectre d'absorption, c'est à dire toutes les valeurs prises par k en fonction de la longueur d'onde, puis de faire des statistiques sur les valeurs de k. Ici, et suivant l'idée de l'échantillonnage des transitions développées notamment dans Galtier et al., 2017 et Nyffenegger-Pere et al., 2024, et rappelées dans le document disponible ici, on ne calcule pas le spectre d'absorption : puisque le coefficient d'absorption à une fréquence donnée est une somme sur les contributions de toutes les raies du spectre, on peut exprimer k comme l'espérance d'une variable aléatoire H qui représente la contribution d'une raie donnée. On peut alors emboiter les espérances (moyenne sur la bande étroite et somme sur les raies) et estimer k_barre ou k^2_barre sans jamais avoir calculer explicitement toutes les valeurs de k_nu. On pourra revenir aux fondamentaux présentés dans le TP Monte Carlo théorique voir ici pour travailler en profondeur cette idée d'emboitement des espérances.

Le programme sln-stat attend en entrée une plage de nombres d'onde, une structure hiérarchique permettant d'accélérer l'échantillonage des raies, et les fichiers HITRAN décrivant les paramètres des raies. Il est possible de réutiliser les fichiers produits au cours du TP1, dans le répertoire TP1/data/ : mixture_tree.bin, molparams.txt et mixture_lines.shtr. Alternativement, il est possible d'en générer de nouveaux en suivant la procédure décrite dans le TP1. Pour utiliser les données du TP1, on va copier le répertoire data/ du TP1 :

cp -r ../TP1/data/ .

Syntaxe :

sln-stat [-dhsv] [-n nrealisations] [-t threads] \
-S nu_min,nu_max -a accel_struct -m molparams \
-l lines

Comme pour le TP1, on va définir la plage spectrale et le nombre de réalisations sur lesquelles les statistiques de k seront calculées :

nu_range=100,125
samples_count=100000

Puis la commande ci-dessous va produire des estimations de E[K_nu] et E[K_nu^2] et les stocker dans les variables bash k et k_2.

read -r k k_2 <<<$(\
sln-stat -S $nu_range -n $samples_count -a data/mixture_tree.bin \
  -m data/molparams.txt -sl data/mixture_lines.shtr \
  | cut -d ":" -f2 | sed 's/+\/-.*$//' | tr -d '\n'\
)

Partie 2 : Estimation des paramètres de Malkmus.

Pour estimer les paramètres de Malkmus k_barre et phi, on va utiliser les résultats de Domoto (1974) et Dufresne et al. (1997) (voir le document sous ${TP_STAR_LINE}/Doc_spectral/dufresne_fournier_grandpeix_inverse_gaussian_kdistributions.pdf également disponible en ligne). Dans son papier de 1974, Domoto exprime la distribution de k des spectres sous-jacents au modèle de Malkmus comme une fonction analytique que Dufresne et al. reconnaissent comme étant la densité de probabilité d'une loi inverse-gaussienne. On peut ainsi faire le lien entre les paramètres de l'inverse-gaussienne, les premiers moments de la distribution et les paramètres du modèle de Malkmus.

Passer au tableau pour retrouver les expressions de l'espérance et variance de K en fonction des paramètres du modèle de Malkmus, puis les paramètres de Malkmus en fonction de E[K_nu] et E[K_nu^2].

On trouve :

E[k] = k_barre
var[k] = k_barre^2/phi

En utilisant l'identité var[k] = E[k^2] - E[k]^2, on obtient :

k_barre = E[K_nu]
phi = E[K_nu]^2/(E[K_nu^2]-E[K_nu]^2)

On peut donc obtenir le paramètre phi de Malkmus à l'aide de la commande suivante :

phi=$(awk -v k=$k -v k_2=$k_2 \
  'BEGIN { print k^2/(k_2-k^2) }')

Enfin, on peut comparer le modèle de Malkmus dont les paramètres ont été estimés aux données générées par sln-slab dans le TP précédent :

ln -sf ../TP1/results.dat .
gnuplot -p -e "k=$k" -e "phi=$phi" \
  gnuplot_scripts/slab_transmissivity.gp