[精讚] [會員登入]
546

[PERL] 常態分配亂數產生函數

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

此文完整連結 http://n.sfs.tw/10536

複製連結 [PERL] 常態分配亂數產生函數@新精讚
(文章歡迎轉載,務必尊重版權註明連結來源)
2017-01-09 11:53:50 最後編修
2017-01-03 01:43:56 By 張○○
 

一般而言我們使用程式內建的 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

你可能感興趣的文章

[PERL] Perl 不立即輸出的列印緩衝區問題 解決Perl 不立即輸出而是最後一次輸出的列印緩衝區問題

[PERL] 使用CPAN安裝模組 在Linux 上,CPAN 可以用來安裝或管理 perl 的模組,此文教你怎麼做。

自行撰寫syslog server建立資訊安全控管中心#4 -- 過濾條件設定 利用PERL將syslog收攏到資料庫的程式,過濾條件設定

[PERL] 18-套件及模組 套件和模組入門

[PERL] 03-條件式判斷 perl的條件式判斷

[PERL] 16-字串取代和置換 Perl 字串比對及置換

[PERL] 簡易檢查網頁記錄檔ip來源統計 利用PERL來檢查網頁記錄檔ip來源統計的簡易程式

Perl 計算經過的時間 Perl 計算程式執行經過的時間

PERL的真值和假值(布林值) 整理Perl中的判斷真假的結果

[PERL] 06-運算子 #2 PERL的運算子介紹,總共有21種

我有話要說


限制:留言最高字數1000字,超過部分會被截掉。請注意:留言不可帶有網址,會被濾掉。 限制:未登入訪客,每則留言間隔需超過10分鐘,每日最多5則留言。

訪客留言

[無留言]

隨機好文

[jQuery] 利用load()來達成ajax的寫法 jQuery中利用load()來達成ajax的寫法,也有人稱他是假的ajax,作法就是..

為什麼要重造輪子? 什麼輪子?造什麼輪子?我為什麼要重造輪子?

TFTP Server 安裝及使用 讓設備的網路設定檔或是韌體經由TFTP拷備出來,操作的方法

NETCRAFT發現你的網站及作業系統 NETCRAFT可以發現你的網站及作業系統

精讚的版面變化 ▓此文僅作為舊文的記錄▓ 這篇文章為了紀念改版完成而撰寫。 原本的部落格是民國97年的作品,那時還是用舊有的技術來寫,很