Estimation du maximum de vraisemblance (MLE) pour la distribution normale avec PROC IML

Niveau de difficulté
Débutant
Publié le :
Le script implémente la fonction de log-vraisemblance pour une distribution normale dans un module SAS©/IML `LogLik`. Il lit la variable `SepalWidth` du dataset `Sashelp.Iris` et utilise la routine d'optimisation non linéaire `nlpnra` pour trouver les estimations du maximum de vraisemblance pour la moyenne (mu) et l'écart-type (sigma). Enfin, il génère un histogramme de la variable et superpose la courbe de densité normale estimée à l'aide de `PROC SGPLOT`, appelée depuis `PROC IML` via des blocs `SUBMIT`/`ENDSUBMIT` pour illustrer l'ajustement du modèle aux données.
Analyse des données

Type : SASHELP


Les données sont tirées du dataset interne `Sashelp.Iris`, spécifiquement la variable `SepalWidth`. Ce jeu de données est disponible par défaut dans l'environnement SAS.

1 Bloc de code
Macro Variables
Explication :
Définit des variables macro pour le nom du jeu de données source (`DSName`) et le nom de la variable d'intérêt (`VarName`). Cela permet une réutilisation facile et une modification centralisée si le jeu de données ou la variable à analyser change.
Copié !
1%let DSName = Sashelp.Iris;
2%let VarName = SepalWidth;
3 
2 Bloc de code
PROC IML
Explication :
Ce bloc `PROC IML` est le cœur de l'analyse. Il définit un module `LogLik` qui calcule la fonction de log-vraisemblance pour une distribution normale, prenant en entrée les paramètres (moyenne `mu` et écart-type `sigma`) et la matrice de données `x`. Les données sont lues dynamiquement depuis `Sashelp.Iris` via les variables macro. Le script initialise les valeurs de départ pour les paramètres (`p`), définit les contraintes (`con`) et les options d'optimisation (`opt`). La fonction `nlpnra` est ensuite appelée pour effectuer l'optimisation non linéaire et trouver les estimations du maximum de vraisemblance (`MLE`), qui sont ensuite affichées.
Copié !
1PROC IML;
2/* write the log-likelihood function for Normal dist */
3start LogLik(param) global (x);
4 mu = param[1];
5 sigma2 = param[2]##2;
6 n = countn(x); /* count nonmissing values */
7 f = -n/2*log(sigma2) - 1/2/sigma2*ssq(x-mu);
8 return ( f );
9finish;
10 
11/* read data */
12use &DSName;
13read all var {&VarName} into x;
14close &DSName;
15/* Specify initial guess for parameters that is close to
16 p = mean(x) || std(x); */
17p = {35 5.5};
18 
19/* mu-sigma constraint matrix */
20con = { . 0, /* lower bounds: -infty < mu; 0 < sigma */
21 . .}; /* upper bounds: mu < infty; sigma < infty */
22opt = {1, /* find maximum of function */
23 0}; /* print iteration history? 0=no; 1=yes */
24call nlpnra(rc, MLE, "LogLik", p, opt, con);
25PRINT MLE[c={"mu" "sigma"}]; /* maximum likelihood estimates */
3 Bloc de code
PROC SGPLOT (via SUBMIT/ENDSUBMIT)
Explication :
Ce bloc démontre comment intégrer d'autres procédures SAS (ici, `PROC SGPLOT`) dans un programme IML en utilisant les instructions `SUBMIT` et `ENDSUBMIT`. Il utilise les estimations `mu` et `sigma` obtenues précédemment pour superposer une courbe de densité normale sur un histogramme de la variable `SepalWidth`, permettant une visualisation de l'ajustement du modèle. L'instruction `quit;` met fin à la procédure IML.
Copié !
1/* Optional: plot histogram with density overlay.
2 The SUBMIT and ENDSUBMIT statements enable you to call
3 SAS procedures from within a SAS/IML program. You can
4 pass parameters to the procedure. See
5 R. Wicklin, "Passing values from PROC IML into SAS procedures",
6 The DO Loop blog, published Jun 3, 2013.
7 http://blogs.sas.com/content/iml/2013/06/03/passing-values-into-procedures/
8*/
9submit mu=(MLE[1]) sigma=(MLE[2]);
10 PROC SGPLOT DATA=&DSName;
11 histogram &VarName;
12 density &VarName / type=normal(mu=&mu sigma=&sigma);
13 RUN;
14endsubmit;
15 
16QUIT;
Ce matériel est fourni "tel quel" par We Are Cas. Il n'y a aucune garantie, expresse ou implicite, quant à la qualité marchande ou à l'adéquation à un usage particulier concernant le matériel ou le code contenu dans les présentes. We Are Cas n'est pas responsable des erreurs dans ce matériel tel qu'il existe maintenant ou existera, et We Are Cas ne fournit pas de support technique pour celui-ci.
Informations de Copyright : Example taken from R. Wicklin, 'Maximum likelihood estimation in SAS/IML', The DO Loop blog, published Oct 12, 2011. URL: http://blogs.sas.com/content/iml/2011/10/12/maximum-likelihood-estimation-in-sasiml/