一般而言我們使用程式內建的 rand 函數取得的亂數是平均分配的，但如果需要常態分配的亂數產生器該怎麼做？

因此我採用Perl改寫亂數產生函式，使其產生1,000,000筆亂數，在我的SERVER上跑的效能極佳，只需要大約4.6秒即可以完成(HP DL380G3 Xeon3.0x2)。

以下是測試的結果，均值為6，標準差為1，我把每隔1統計一筆資料：
以下是程式產生亂數後再分群統計數量，由結果看的確是符合常態分配的函數：

Program Start....
1 3
2 235
3 5919
4 60437
5 242095
6 382825
7 241478
8 60621
9 6177
10 209
11 1
Elapsed time=00:00:04.671079

至於實測10,000,000筆的資料，在 linux OS和INTEL雙CPU(非雙核喔)表現上(沒有寫多執行緒，也許沒法最佳化)，所費的時間毫無意外是1,000,000筆的10倍。這裡如果拿 IA64 來測的話，不知效果如何… 。

Program Start....
1 31
2 2231
3 59565
4 605552
5 2417906
6 3831267  <==均值
7 2415974
8 605805
9 59397
10 2245
11 27
Elapsed time=00:00:46.395085[

以下是 Perl的程式碼片段，分群組只花了三行：

#!/usr/bin/perl -w

use Time::Elapse;
print "Program Start....\n";
$avg = 6; # 均值
$sd =1; # 標準差
$num = 10000000; # 產生組數
$dif = 1; # 統計柱狀圖區間幾個 $sd
# ====== 以上是參數的設定部分

my %h=(); # 分組的HASH
Time::Elapse->lapse(my $now);
for(1..$num){
    do{
        $v1 = (2 * rand 1) - 1;
        $v2 = (2 * rand 1) - 1;
        $r = $v1 **2 + $v2 **2;
    }while($r >= 1);
    $fac = sqrt(-2 * log($r)/ log(exp(1)) / $r);
    $gauss = $v2 * $fac;
    $vr= $gauss * $sd + $avg;     # $vr 就是得到的亂數
    # 計算群組中間值
    $sddif = $sd*$dif;
    $hid= int( ( $vr + $sddif /2 ) / $sddif ) * $sddif;
    # 置入hash
    $h{$hid}++;
}

foreach $key (sort {$a<=>$b} keys %h) {
    print "$key \t\t$h{$key}\n";
}
print "Elapsed time=". $now;

如果只要一個「常態分配」的亂數產生函數而不要像上面的統計分析，我把他抽出來變副程式：

#!/usr/bin/perl -w
use strict;

#&mean_dis_fun(均值, 標準差);
print &mean_dis_fun(5, 1);

sub mean_dis_fun{
    my $avg =shift;
    my $sd = shift;
    my ($r, $v1, $v2);
    do{
        $v1 = (2 * rand 1) - 1;
        $v2 = (2 * rand 1) - 1;
        $r = $v1 **2 + $v2 **2;
    }while($r >= 1);
    my $fac = sqrt(-2 * log($r)/ log(exp(1)) / $r);
    my $gauss = $v2 * $fac;
    return $gauss * $sd + $avg;
}

原文撰於 2008-10-11，重新整理於 2009-11-15