使用Fortran在一個範圍內生成正態分佈的數據

Question

我想根據給定的均值和給定的上限和下限生成一個正態分佈數據的數組。

我發現使用Marsaglia極座標法（a.k.a polar form of the Box-Müller transformation）的函數here。

我該如何修改它以在均值附近的-10和+10整數值內生成點？

function normal(mean, sigma) 

    ! mean : mean of distribution 
    ! sigma : number of standard deviations 

    implicit none 

    integer, parameter:: b8 = selected_real_kind(14) 
    real(b8), parameter :: pi = 3.141592653589793239_b8 
    real(b8) normal 
    real rand_num 
    real(b8) tmp 
    real(b8) mean 
    real(b8) sigma 
    real(b8) fac 
    real(b8) gsave 
    real(b8) rsq 
    real(b8) r1 
    real(b8) r2 
    integer flag 
    save flag 
    save gsave 
    data flag /0/ 

    if (flag.eq.0) then 
    rsq=2.0_b8 

    do while(rsq.ge.1.0_b8.or.rsq.eq.0.0_b8) ! new from for do 
     call random_number(rand_num) 
     r1=2.0_b8*rand_num-1.0_b8 
     call random_number(rand_num) 
     r2=2.0_b8*rand_num-1.0_b8 
     rsq=r1*r1+r2*r2 
    enddo 

    fac=sqrt(-2.0_b8*log(rsq)/rsq) 
    gsave=r1*fac 
    tmp=r2*fac 
    flag=1 
    else 
    tmp=gsave 
    flag=0 
    endif 

    normal=tmp*sigma+mean 

    return 
    endfunction normal 

Answer 1

從截斷的正態分佈中採樣有多種不同的方法。 「最好」的方式取決於你的均值和方差。

一個簡單的方法就是plain rejection：產生偏差，如果它太大或太小，則將其丟棄並重復。如果你的參數是這樣的，你不拒絕許多這是一個好方法。

對於指定的程序，做這種拒絕很簡單：堅持整個計算部分（包括flag測試）爲do ... end do環路以if (abs(normal).lt.10) exit應該這樣做。

[這不是優雅的，可能需要額外的簿記，樣本是以成對的方式生成的。但是，如果拒絕發生很少，那也不算太壞。整理是作爲一個練習讀者。]

這就是如何可以改變例程，但作爲@george評論，這可能不是最好的方法，即使拒絕。當然，有一個叫做normal的函數，它需要mean和sigma，並返回一個不正常分佈的偏差，這可能會令人困惑。

function truncated_normal(mu, sigma, lim_a, lim_b) !Very side-effecty 
    ... declarations ... 
    if (... test for wanting the rejection method ...) then 
    do 
     truncated_normal = normal(mu, sigma) 
     if (truncated_normal.gt.lim_a.and.truncated_normal.lt.lim_b) exit 
    end do 
    else 
    ... the other magical method when rejection is too likely ... 
    end if 
end