logo标签列表

在Perl中，如何计算正态分布下某一点的概率？

perlrstatisticsprobability
3

Perl中是否有一个包可以计算每个给定点的概率分布高度。例如，在R中可以这样做：

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

点x=0的概率符合均值为4，标准差为10的正态分布，为0.0368。我查看了Statistics::Distribution，但它并没有提供这个函数。
2
正态分布中任意一点的概率当然是零。你想要计算什么？ - Eduardo Leoni
@EL：我不是指“任意/随机”的点，而是一个“特定”的点。 - neversaint
3
正态分布是连续的，因此任何点（无论是否给定）的概率都是零。也许你想要密度？这就是dnorm中的“d”代表的意思。 - Eduardo Leoni
7个回答
8

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（Debian软件包r-mathlib提供）编写Perl扩展来获取这些函数？这不需要R解释器。

否则，您可以尝试GNU GSL或Cephes库以访问这些特殊函数。

CPAN上已经有一个可以使用R的模块。虽然它很混乱，但我过去曾经让它工作过：http://search.cpan.org/~gmpassos/Statistics-R-0.02/ - tsee
2
统计学中的分布函数（如 pnorm）在 Statistics::Distributions 中被称为 uprob。1-uprob（（0-4）/10）应该给你 ~ 0.34（我没有安装它来确认这一点）。不过，我没有密度函数。 - Eduardo Leoni
4
为什么不使用类似以下方式（我用R语言编写，但也可以使用Statistics::Distribution模块的perl语言实现）：
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] 修复了复制粘贴错误。

@EL：谢谢。当最终结果为负数时，您将如何调整您的方法？例如，x=2.02，平均值=2，标准差=0.24。您的方法会得出-2.880624e-05。 - neversaint
在您的最后一个例子中，我的机器给出了不同的结果：dn(2.02,2,.24); [1] 1.656469。我正在使用R版本2.9.2。 - neversaint
1
@foolishbrat：那是我的错误。1.65才是正确答案（与dnorm的答案一致）。对于造成的困惑，非常抱歉。 - Eduardo Leoni
人们通常在dnorm大于1时会做什么，就像这个例子一样？ - neversaint
2
@foolishbrat：我认为你再次混淆了概率（其范围在0和1之间）和概率密度（其范围不在此之间）。正如其他人指出的那样，你可能想要累积分布函数；但由于你没有告诉我们你想做什么，所以我们无法知道。你还应该查阅一本统计学入门书籍。 - Eduardo Leoni
3

如果您真的需要密度函数，为什么不直接使用它：

$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

我认为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
使用 Perl 的 Statistics::Distributions 模块，您可以通过以下方式实现：
#!/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”的结果。

1
正如其他人指出的那样，您可能需要累积分布函数。可以通过误差函数（平均值平移并缩放为正态分布的标准差）获得此函数，该函数存在于标准数学库中，并且通过Math::Libm在Perl中进行访问。
0

以下是如何使用 CPAN 中的 Math::SymbolicX::Statistics::Distributions 模块在 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

这应该从Dirk的答案中给你 ~0.344578258389676。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接