如何在Perl中计算给定正态分布的点的概率?

How can I compute the probability at a point given a normal distribution in Perl?

提问人:neversaint 提问时间:9/4/2009 最后编辑:brian d foyneversaint 更新时间:9/28/2016 访问量:9325

问:

Perl 中是否有一个包可以让你计算每个给定点的概率分布高度。例如,这可以在 R 中以这种方式完成:

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

也就是说,点 x=0 落入正态分布的概率为 0.0368,均值 = 4,sd=10。 我看了一下 Statistics::D istribution,但它并没有给出很多 函数来做到这一点。

Perl R 统计 概率

评论

2赞 Eduardo Leoni 9/4/2009
正态分布中任何点的概率当然为零。你想计算什么?
0赞 neversaint 9/4/2009
@EL:我指的不是“任何/随机”,而是“给定”的观点。
3赞 Eduardo Leoni 9/4/2009
法线是连续的,因此任何点(给定或不给定)的概率为零。也许你想要密度?(这就是 dnorm 中的“d”代表。

答:

8赞 Dirk is no longer here 9/4/2009 #1

dnorm(0, mean=4, sd=10) 不会给出此类点发生的概率。引用维基百科关于概率密度函数的话

在概率论中,概率 密度函数 (PDF) - 通常引用 作为概率分布 函数1 - 或随机密度 变量是一个函数,用于描述 每个概率的密度 点在样本空间中。这 随机变量的概率 落在给定集合内的下式由 其密度在 设置。

你提到的概率是

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

或从 N(4, 10) 分布中获得等于或小于 0 的值的几率为 34.46%。

至于你的Perl问题:如果你知道如何在R中做到这一点,但需要从Perl中得到它,也许你需要编写一个基于R的libRmath的Perl扩展(在Debian中由r-mathlib软件包提供)来将这些函数带到Perl?这不需要 R 解释器。

否则,您可以尝试使用 GNU GSL 或 Cephes 库来访问这些特殊功能。

评论

0赞 tsee 9/4/2009
CPAN 上已经有一个可以使用 R 的模块。这是一团糟,但我过去可以让它工作:search.cpan.org/~gmpassos/Statistics-R-0.02
2赞 Eduardo Leoni 9/4/2009
Statistics::D istributions 中的分布函数(如 pnorm)是 uprob。1-uprob((0-4)/10) 应该给你 ~ 0.34(我没有安装它来确认这一点。不过,我没有密度函数。
4赞 Eduardo Leoni 9/4/2009 #2

为什么不按照这些思路(我是用 R 写的,但它可以在 perl 中使用 Statistics::D istribution 完成):

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

[编辑:1]我应该提一下,这种近似值在远尾可能会变得非常可怕。这可能重要,也可能无关紧要,具体取决于您的应用程序。

[编辑:2] @foolishbrat 将代码转换为函数。结果应该始终是积极的。也许你忘记了在perl模块中你提到函数返回上限概率1-F,而R返回F?

[编辑:3] 修复了复制和粘贴错误。

评论

0赞 neversaint 9/4/2009
@EL:谢谢。当最终结果是负面的时,你会如何调整你的方法。例如,x=2.02,平均值=2,sd=0.24。您的方法将给出 -2.880624e-05。
0赞 neversaint 9/7/2009
@EL:在上一个示例中。我的机器给出了不同的结果:dn(2.02,2,.24);[1] 1.656469.我使用的是 R 版本 2.9.2。
1赞 Eduardo Leoni 9/7/2009
@foolishbrat:这是我的错误。~1.65 是正确的。(并同意 Dnorm 的回答。很抱歉造成混乱。
0赞 neversaint 9/7/2009
@EL:当 dnorm 大于 1 时,人们通常会怎么做?
2赞 Eduardo Leoni 9/7/2009
@foolishbrat:再说一次,我认为你混淆了概率(介于 0 和 1 之间)和概率密度(不是)。就像其他人指出的那样,您可能想要累积分布函数;但是我们没有办法知道,因为你没有告诉我们你想做什么。您还应该查阅统计书的介绍。
0赞 tsee 9/4/2009 #3

以下是如何使用 CPAN 的 Math::SymbolicX::Statistics::D istributions 模块在 Perl 中使用 R 执行相同的操作:

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

该模块中的 normal_distribution() 函数是函数的生成器。$norm将是可以修改的 Math::Symbolic (::Operator) 对象。例如,使用 implement,在上面的示例中,它将两个参数变量替换为常量。

然而,请注意,正如 Dirk 所指出的,你可能想要正态分布的累积函数。或者更一般地说,是一定范围内的积分。

不幸的是,Math::Symbolic 不能以符号方式进行积分。因此,您必须求助于Math::Integral::Romberg之类的数值积分。(或者,在 CPAN 中搜索错误函数的实现。这可能很慢,但仍然很容易做到。将以下内容添加到上面的代码片段中:

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

这应该给你 ~0.344578258389676 来自 Dirk 的答案。

1赞 Jouni K. Seppänen 9/5/2009 #4

正如其他人所指出的,您可能想要累积分布函数。这可以通过误差函数(按均值平移,按正态分布的标准差缩放)获得,该函数存在于标准数学库中,并且可以通过 Math::Libm 在 Perl 中访问。

3赞 Eonwe 10/25/2010 #5

如果你真的想要密度函数,为什么不直接使用它:

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

它给出的 1.65649768474891 与 R 中的 dnorm 大致相同。

2赞 Jonathan Ledlie 3/20/2012 #6

我不认为 Jouni 是完全正确的。这似乎给出了一个合理的 PDF 版本(如果您只想要一个特定的 x-y 点,请提取循环的中间部分):

!/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);
}
1赞 Kit 8/27/2013 #7

使用 Perl 的 Statistics::D istributions,您可以通过以下方式实现此目的:

#!/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";

“点 0 处的概率分布高度 = 0.34458”的结果