Comment puis-je calculer la probabilité à un point donné une distribution normale en Perl?

Question

Y at-il un paquet en Perl qui vous permet de calculer la hauteur de la distribution de probabilité à chaque point donné. Par exemple, cela peut se faire en R cette façon:

> dnorm(0, mean=4,sd=10) 
> 0.03682701 

A savoir la probabilité du point x = 0 tombe dans une distribution normale, avec une moyenne = 4 et sd = 10, est 0,0368. J'ai regardé Statistics::Distribution mais cela ne donne pas cette fonction très pour le faire.

Answer 1

Pourquoi ne pas quelque chose dans ce sens (je vous écris en R, mais il pourrait être fait en perl avec Statistique :: Distribution):

dn <- function(x=0 # value 
       ,mean=0 # mean 
       ,sd=1 # sd 
       ,sc=10000 ## scale the precision 
       ) { 
    res <- (pnorm(x+1/sc, mean=mean, sd=sd)-pnorm(x, mean=mean, sd=sd))*sc 
    res 
} 
> dn(0,4,10,10000) 
0.03682709 
> dn(2.02,2,.24) 
1.656498 

[edit: 1] Je dois mentionner que cette approximation peut devenir assez horrible à l'extrême. cela pourrait ou non importer en fonction de votre application.

[edit: 2] @foolishbrat Transformé le code en fonction. Le résultat devrait toujours être positif. Vous oubliez peut-être que dans le module perl vous mentionnez que la fonction renvoie la probabilité supérieure 1-F et que R renvoie F? [Edit: 3] Correction d'une erreur de copier-coller.

[edit: 3]

Answer 2

dnorm (0, moyenne = 4, sd = 10) ne pas vous donner la probabilité thr d'un tel point se produisant. Pour citer Wikipedia sur probability density function

En théorie des probabilités, une probabilité fonction de densité (pdf) -souvent appelé comme une distribution de probabilité fonction 1 densité -ou, d'une variable aléatoire est une fonction qui décrit la densité de probabilité à chaque point dans l'espace d'échantillonnage. La probabilité d'une variable aléatoire appartenant à un ensemble donné est donnée par l'intégrale de sa densité par rapport à l'ensemble .

et la probabilité vous mentionner est

R> pnorm(0, 4, 10) 
[1] 0.3446 

ou une chance 34,46% d'obtenir une valeur égale ou inférieure à 0 d'une distribution N (4, 10). En ce qui concerne votre question Perl: Si vous savez comment faire en R, mais que vous avez besoin de Perl, vous devrez peut-être écrire une extension Perl basée sur R libRmath (fournie par Debian par le paquetage r-mathlib) pour obtenir ces fonctions à Perl? Cela ne nécessite pas l'interpréteur R. Dans le cas contraire, vous pouvez essayer les bibliothèques GNU GSL ou les céphes pour accéder à ces fonctions spéciales.

Answer 3

Voilà comment vous pouvez faire la même chose que vous faites avec R en Perl utilisant le module Math::SymbolicX::Statistics::Distributions de CPAN:

use strict; use warnings; 

use Math::SymbolicX::Statistics::Distributions qw/normal_distribution/; 

my $norm = normal_distribution(qw/mean sd/); 
print $norm->value(mean => 4, sd => 10, x => 0), "\n"; 

# curry it with the parameter values 
$norm->implement(mean => 4, sd => 10); 
print $norm->value(x => 0),"\n"; # prints the same as above 

La fonction normal_distribution() à partir de ce module est un générateur de fonctions. $ norm sera un objet Math::Symbolic (:: Operator) que vous pouvez modifier.Par exemple avec mettre en œuvre, qui, dans l'exemple ci-dessus, remplace les deux variables de paramètre par des constantes. Notez cependant, comme l'a souligné Dirk, que vous voulez probablement la fonction cumulative de la distribution normale. Ou plus généralement l'intégrale dans une certaine gamme.

Malheureusement, Math :: Symbolic ne peut pas effectuer d'intégration symboliquement. Par conséquent, vous devrez recourir à l'intégration numérique avec les goûts de Math::Integral::Romberg. (Vous pouvez également rechercher dans le CPAN une implémentation de la fonction d'erreur.) Cela peut être lent, mais c'est toujours facile à faire. Ajoutez ceci à l'extrait ci-dessus:

use Math::Integral::Romberg 'integral'; 

my ($int_sub) = $norm->to_sub(); # compile to a faster Perl sub 
print $int_sub->(0),"\n"; # same number as above 

print "p=" . integral($int_sub, -100., 0) . "\n"; 
# -100 is an arbitrary, small number 

Cela devrait vous donner la ~ ,344578258389676 de la réponse de Dirk.

Answer 4

Comme d'autres l'ont souligné, vous voulez probablement la fonction de distribution cumulative. Ceci peut être obtenu via le error function (décalé par la moyenne et mis à l'échelle par l'écart-type de votre distribution normale), qui existe dans la bibliothèque mathématique standard et est rendu accessible en Perl par Math::Libm.

Answer 5

Si vous voulez vraiment la fonction de densité, pourquoi ne pas utiliser directement:

$pi = 3.141593; 
$x = 2.02; 
$mean = 2; 
$sd = .24; 
print 1/($sd * sqrt(2*$pi)) * exp(-($x-$mean)**2/(2 * $sd**2)); 

Il donne 1,65649768474891 environ le même que dnorm R.

Answer 6

Je ne pense pas Jouni est tout à fait raison. Cela semble donner une version raisonnable du PDF (extrait au milieu de la boucle si vous voulez juste un point xy spécifique):

!/usr/bin/perl 

use strict; 
use Getopt::Std; 
use POSIX qw(ceil floor); 

# Usage 
# Outputs normal density function given a mean and sd 
# -s standard deviation 
# -m mean 
# -n normalization factor (multiply result by this amount), optional 

my %para =(); 
getopts('s:m:n:', \%para); 
if (!exists ($para{'s'}) || !exists ($para{'m'})) { 
    die ("mean and standard deviation required"); 
} 

my $norm = 1.0; 
if (exists ($para{'n'})) { 
    $norm = $para{'n'}; 
} 

my $sd = $para{'s'}; 
my $mean = $para{'m'}; 

my $start = floor($mean - ($sd * 5)); 
my $end = ceil($mean + ($sd * 5)); 

my $pi = 3.141593; 

my $var = $sd**2; 

for (my $x = $start; $x < $end; $x+=0.1) { 
    my $e = exp(-1 * (($x-$mean)**2)/(2*$var)); 
    my $d = sqrt($var) * sqrt(2*$pi); 
    my $y = 1.0/$d*$e * $norm; 
    printf ("%5.5f %5.5f\n", $x, $y); 
} 

Answer 7

En utilisant les statistiques de Perl :: Les distributions, vous pouvez y parvenir avec:

#!/usr/bin/perl 

use strict; use warnings; 
use Statistics::Distributions qw(uprob); 

my $x  = 0; 
my $mean = 4; 
my $stdev = 10; 

print "Height of probablility distribution at point $x = " 
    . (1-uprob(($x-$mean)/$stdev))."\n"; 

Résultats avec « Hauteur de la distribution de probablility au point 0 = 0,34458 »