Le laboratoire de science des données
Les fonctions LogBeta et LogGamma utilisant C#
La bibliothèque .NET (anciennement appelée .NET Core) n’a pas de fonctions intégrées pour les analyses statistiques classiques. Mais il est possible d’implémenter de telles fonctions à partir de zéro. Cet article présente des implémentations C# de trois des fonctions statistiques classiques les plus importantes : la fonction log-bêta, la fonction log-gamma et la fonction bêta incomplète régularisée.
La fonction bêta incomplète régularisée s’écrit généralement Ix(a,b) ou I(x; a,b). Il est utilisé pour de nombreux problèmes, y compris l’analyse statistique de la variance (ANOVA). La fonction Ix(a,b) appelle la fonction log-beta. La fonction log-beta appelle la fonction log-gamma.
Une bonne façon de voir où va cet article est de jeter un coup d’œil à la capture d’écran d’un programme de démonstration dans Figure 1. La première partie de la sortie montre la valeur de I(x; a, b) pour x = 0,5, a = 4,0, b = 7,0, soit 0,828125. La fonction bêta incomplète régularisée n’a pas d’interprétation évidente. Il vaut mieux la considérer comme une fonction mathématique abstraite qui accepte une valeur x comprise entre 0,0 et 1,0, et des valeurs a et b positives. La valeur de retour est un nombre compris entre 0,0 et 1,0, c’est pourquoi la fonction est appelée “régularisée”.
La deuxième partie de la sortie dans Figure 1 affiche la valeur de LogBeta(3.0, 5.0), qui est -4,653960. Encore une fois, la fonction n’a pas d’interprétation évidente et il vaut mieux la considérer simplement comme une fonction qui accepte deux valeurs positives. Il existe une simple fonction Beta(a,b) mais la valeur de retour de Beta() peut être très grande et il est donc courant d’utiliser LogBeta() qui est littéralement le journal de la fonction Beta().
La troisième partie de la sortie dans Figure 1 affiche la valeur de LogGamma(4.0), qui est de 6,0000. La fonction simple Gamma() est une généralisation de la fonction Factorial(). Si n est un entier, alors Gamma(n) = Factoriel(n-1). Cependant, Gamma(z) est défini pour un z réel positif. Étant donné que la valeur de retour de la fonction Gamma(z) peut être astronomique, il est courant d’utiliser LogGamma(z).
Cet article suppose que vous avez des compétences de programmation intermédiaires ou supérieures avec le langage C #, mais ne suppose pas que vous savez quoi que ce soit sur les fonctions bêta, gamma ou bêta incomplètes régularisées. Le code de démonstration complet est présenté dans cet article et est également disponible dans le téléchargement qui l’accompagne. Toutes les vérifications d’erreur normales ont été supprimées pour garder les idées principales aussi claires que possible.
Pour créer le programme de démonstration, j’ai utilisé Visual Studio 2022 (l’édition communautaire gratuite) avec .NET 5 (pour éviter le nouveau modèle de console .NET 6 ennuyeux). Cependant, la démo n’a pas de dépendances significatives, donc toutes les versions relativement récentes de VS et .NET peuvent être utilisées.
La fonction LogGamma()
Le graphique de la fonction simple Gamma() est affiché dans Figure 2. Notez que Gamma(n) = Factorial(n-1), par exemple, Gamma(4.0) = Factorial(3) = 3 * 2 * 1 = 6.0.
La fonction LogGamma(z) ne peut pas être calculée exactement lorsque z n’est pas un entier et il existe donc plusieurs algorithmes d’approximation. Le programme de démonstration utilise ce qu’on appelle l’algorithme Lanczos g=5, n=7. La mise en œuvre est illustrée dans Liste 1.
Liste 1 :
La fonction LogGamma()
static double LogGamma(double z)
{
// Lanczos approximation g=5, n=7
double[] coef = new double[7] { 1.000000000190015,
76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2, -0.5395239384953e-5 };
double LogSqrtTwoPi = 0.91893853320467274178;
if (z < 0.5)
return Math.Log(Math.PI / Math.Sin(Math.PI * z)) -
LogGamma(1.0 - z);
double zz = z - 1.0;
double b = zz + 5.5; // g + 0.5
double sum = coef[0];
for (int i = 1; i < coef.Length; ++i)
sum += coef[i] / (zz + i);
return (LogSqrtTwoPi + Math.Log(sum) - b) +
(Math.Log(b) * (zz + 0.5));
}
Le programme de démonstration appelle la fonction LogGamma() en utilisant ces instructions :
double z = 4.0;
double lg = LogGamma(z);
Console.WriteLine("LogGamma(4.0) = ");
Console.WriteLine(lg.ToString("F6"));
double g = Math.Exp(lg);
Console.WriteLine("Gamma(4.0) = ");
Console.WriteLine(g.ToString("F6"));
Le calcul de Gamma(4.0) en appelant Math.Exp() sur LogGamma(4.0) est fait juste pour montrer que la fonction LogGamma() n'est que le logarithme de la fonction gamma.
La démo ne fait aucune vérification d'erreur. Dans un scénario non démo, vous souhaiterez probablement vérifier la valeur du paramètre d'entrée pour vous assurer qu'elle est positive.
La fonction LogBeta()
La fonction bêta incomplète régularisée I(x; a, b) appelle la fonction LogBeta() qui à son tour appelle la fonction LogGamma(). L'implémentation de démonstration de la fonction LogBeta() est :
static double LogBeta(double a, double b)
{
double logbeta = LogGamma(a) + LogGamma(b) -
LogGamma(a + b);
return logbeta;
}
Le programme de démonstration appelle la fonction LogBeta() comme suit :
double lb = LogBeta(3.0, 5.0);
Console.WriteLine("LogBeta(3.0, 5.0) = ");
Console.WriteLine(lb.ToString("F6"));
Mathématiquement, beta(a, b) = (gamma(a) * gamma(b)) / gamma(a + b). Si vous prenez le log() des deux côtés de l'équation et utilisez les faits que log(x / y) = log(x) - log(y) et log(x * y) = log(x) + log( y) vous obtenez le résultat dans l'implémentation de la fonction LogBeta().
La fonction bêta incomplète régularisée
La fonction bêta incomplète régularisée I(x; a, b) est extrêmement délicate à mettre en œuvre. Le programme de démonstration utilise ce qu'on appelle l'approximation des fractions continues (CF). Les équations d'approximation CF sont présentées dans figure 3.
L'approximation I(x; a, b) CF est le produit mathématique de deux termes principaux. Le terme de gauche implique Beta(a, b). Le terme de droite a un nombre infini de d[i] termes où la valeur de chaque terme dépend du fait que [i] est pair ou impair. La démo implémente la fonction I(x; a, b) comme indiqué dans Liste 2.
Liste 2 :
Fonction bêta incomplète régularisée
static double RegIncBeta(double x, double a, double b)
{
// pick the form of RegIncompleteBeta() that converges best
if (x < (a + 1.0) / (a + b + 2.0))
return RegIncompleteBeta(x, a, b);
else
return 1.0 - RegIncompleteBeta(1 - x, b, a);
}
static double RegIncompleteBeta(double x, double a, double b)
{
// regularized incomplete beta
// calls LogBeta() and helper ContFraction()
double cf = ContFraction(x, a, b);
double logTop = (a * Math.Log(x)) +
(b * Math.Log(1 - x));
double logBot = Math.Log(a) + LogBeta(a, b);
double logLeft = logTop - logBot;
double left = Math.Exp(logLeft);
return left * cf; // should be [0.0, 1.0]
}
static double ContFraction(double x, double a, double b,
int maxTerms = 100)
{
// 1. pre-compute 100 d values
double[] d = new double[maxTerms]; // d[0] not used
int end = maxTerms / 2; // 50
for (int m = 0; m < end; ++m) // [0,49]
{
int i = 2 * m; // even di
int j = i + 1; // odd di
d[i] = (m * (b - m) * x) / ((a + 2 * m - 1) *
(a + 2 * m));
d[j] = -1 * ((a + m) * (a + b + m) * x) / ((a + 2 * m) *
(a + 2 * m + 1));
}
// 2. work backwards
double[] t = new double[maxTerms]; // t[0] not used
t[maxTerms - 1] = 1 + d[maxTerms - 1];
for (int j = maxTerms - 2; j >= 1; --j)
t[j] = 1 + d[j] / t[j + 1];
// ShowVector
// Console.ReadLine();
double cf = 1 / t[1];
return cf;
}
Il existe une fonction RegIncBeta() de niveau supérieur. Il appelle la fonction RegIncompleteBeta() qui fait tout le travail. L'approche wrapper est basée sur la propriété mathématique que I(x; a, b) = 1 - I(x-1; b, a) et elle appelle la version qui converge plus rapidement. La fonction RegIncompleteBeta() appelle une fonction auxiliaire ContFraction() qui calcule le terme CF. Pour éviter de calculer Beta(a, b), qui peut facilement déborder, la fonction RegIncompleteBeta() calcule LogBeta(a, b) puis annule l'opération de journalisation par la fonction d'appel Math.Exp().
Le programme de démonstration appelle la fonction I(x; a, b) en utilisant ces instructions :
double x = 0.5;
double a = 4.0;
double b = 7.0;
double rib = RegIncBeta(x, a, b);
Console.WriteLine("Ix(a,b) for (x=0.5, a=4.0, b=7.0) = ");
Console.WriteLine(rib.ToString("F6"));
La fonction d'assistance ContFraction() pré-calcule 100 d[i] valeurs puis calcule le terme CF en partant du dernier d[i] valeur vers le premier d[i] évaluer. Cette approche peut échouer avec un sous-dépassement ou un débordement numérique pour certaines combinaisons de x, a et b, donc dans un scénario non démo, vous souhaiterez probablement envelopper le code dans un bloc try-catch.
Au lieu d'utiliser 100 d pré-calculés[i] termes pour la fonction d'assistance ContFraction(), vous pouvez examiner la valeur de la fonction RegIncBeta(), en commençant par maxTerms = 2 et en augmentant maxTerms jusqu'à ce que la fonction RegIncBeta() converge vers une valeur stable. Par exemple, avec x = 0,5, a = 4,0, b = 7,0 la fonction RegIncBeta() se stabilise à maxTerms = 10 :
maxTerms RegIncBeta()
=======================
2 0.812500
4 0.828613
6 0.828119
8 0.828125
10 0.828125
Une approche alternative consiste à utiliser ce qu'on appelle l'algorithme de Lentz pour estimer la valeur du terme CF. Au lieu de pré-calculer un nombre arbitraire de d[i] termes, puis en revenant à partir du dernier d[i] valeur, l'algorithme de Lentz utilise quelques astuces mathématiques astucieuses pour aller de l'avant, en vérifiant à chaque itération pour voir si la valeur CF a convergé vers une valeur stable. Cependant, l'algorithme de Lentz est souvent très fragile et pour certains problèmes ne converge pas.
Emballer
Il existe de nombreuses bibliothèques de codes externes qui ont des implémentations de fonctions statistiques complexes telles que bêta, log-bêta, gamma, log-gamma et bêta incomplet régularisé. Par exemple, la bibliothèque de langage Scipy Python contient ces fonctions. Dans de nombreuses situations, vous pouvez utiliser ces bibliothèques externes. Cependant, l'implémentation de fonctions statistiques à partir de zéro à l'aide de C # est utile lorsque vous souhaitez éviter une dépendance externe pour des raisons techniques ou de droits d'auteur.
A propos de l'auteur
Dr. James McCaffrey travaille pour Microsoft Research à Redmond, Wash. Il a travaillé sur plusieurs produits Microsoft, dont Azure et Bing. Jacques est joignable au [email protected].
.