ヒューリンクステクニカルサポート
モンテカルロ法の例
SYSTAT 製品ページ

TOP > モンテカルロ> 例

棄却サンプリング

α が整数のとき、独立した α exponential(β) を加算して gamma(α,β) から無作為標本を生成することができます。 α が整数でない場合は、この簡単な方法を使用することはできません。 SYSTAT の単変量連続型ランダム標本抽出プロシージャを使用して gamma(α,β) から無作為標本を生成することができますが、この例では、別の方法を示します。uniform(0,15)、exponential(2.43)、および gamma ([2.43],2.43/[2.43]) 分布をプロポーザルとみなすことにより棄却サンプリングを使用します。 [2.43] は 2.43 の整数部を表します。

  • 一様密度関数をプロポーザル関数として使用して から無作為標本を生成し、この標本から基本統計量を計算し、モンテカルロ積分法を使用して E(X2) の近似を行います。 次のように入力します。
    IIDMC
    SAVE REJFLAT.SYD
    RJS TARGET='(1/(EXP(LGM(2.43))))*EXP(-X)*X^1.43',
    CONSTANT=4.7250 /SIZE=100000 NSAMPLE=1,
    RSEED=3245425
    PROPOSAL U(0,15)
    GENERATE
    INTEG FUN='X^2'/MC
    ESTIMATE
    USE REJFLAT.SYD
    STATS
    CBSTAT S1/ MAXIMUM MEAN MINIMUM SD VARIANCE N
    DENSITY S1 / HIST

    出力は次のようになります。

    Classical Monte-Carlo Integration estimates for S1 :  
    Empirical Mean : 8.308.
    Standard Error : 0.096.
      S1
    N of cases 100000
    Minimum 0.020
    Maximum 14.395
    Mean 2.426
    Standard Dev 1.557
    Variance 2.424


  • 指数密度関数をプロポーザル関数として使用して から無作為標本を生成し、この標本から基本統計量を計算し、モンテカルロ積分法を使用して E(X2) の近似を行います。 次のように入力します。
    IIDMC
    SAVE REJEXPO.SYD
    RJS TARGET='(1/(EXP(LGM(2.43))))*EXP(-X)*X^1.43',
    CONSTANT= 1.6338/ SIZE=100000 NSAMPLE=1,
    RSEED=534652
    PROPOSAL E(0,2.43)
    GENERATE
    INTEG FUN='X^2'/MC
    ESTIMATE
    USE REJEXPO.SYD
    STATS
    CBSTAT S1/ MAXIMUM MEAN MINIMUM SD VARIANCE N
    DENSITY S1 / HIST

    出力は次のようになります。

    Classical Monte-Carlo Integration estimates for S1 :  
    Empirical Mean : 8.323.
    Standard Error : 0.096.
      S1
    N of cases 100000
    Minimum 0.005
    Maximum 15.603
    Mean 2.429
    Standard Dev 1.556
    Variance 2.422


  • ガンマ密度関数をプロポーザル関数として使用して から無作為標本を生成し、この標本から基本統計量を計算し、モンテカルロ積分法を使用して E(X2) の近似を行います。次のように入力します。
    IIDMC
    SAVE REJGAM.SYD
    RJS TARGET='(1/(EXP(LGM(2.43))))*EXP(-X)*X^1.43',
    CONSTANT=1.1102/SIZE=100000 NSAMPLE=1,
    RSEED=236837468
    PROPOSAL G(2,1.2150)
    GENERATE
    INTEG FUN='X^2'/MC
    ESTIMATE
    USE REJGAM.SYD
    STATS
    CBSTAT S1/ MAXIMUM MEAN MINIMUM SD VARIANCE N
    DENSITY S1 / HIST

    出力は次のようになります。

    Classical Monte-Carlo Integration estimates for S1 :  
    Empirical Mean : 8.339.
    Standard Error : 0.096.
      S1
    N of cases 100000
    Minimum 0.008
    Maximum 15.376
    Mean 2.430
    Standard Dev 1.560
    Variance 2.434


  • ガンマ、指数、および一様分布をプロポーザル関数として使用し、密度関数 から無作為標本を生成します。 ただし、目的の関数から標本を受け入れる確率 1/M はプロポーザル関数によって異なります。

    プロポーザル関数 受け入れ確率
    一様 0.21164021
    指数 0.61207002
    ガンマ 0.90073861

    目的の変量としてプロポーザル標本を受け入れる確率は、プロポーザルの積と定数がターゲット関数にどれだけ近いかに依存します。 ターゲット関数とプロポーザルの積および定数をプロットして、これを観察します。 次のように入力します。
    BEGIN
    FPLOT Y= GDF(X,2.43,1); XMIN=0 XMAX=10 YMIN=0 YMAX=0.4,
    HEIGHT=2 WIDTH=2
    FPLOT Y=4.7250*UDF(X,0,15);XMIN=0 XMAX=10 YMIN=0 YMAX=0.4,
    COLOR=1 HEIGHT=2 WIDTH=2
    END
    BEGIN
    FPLOT Y= GDF(X,2.43,1); XMIN=0 XMAX=10 YMIN=0 YMAX=0.6,
    HEIGHT=2 WIDTH=2
    FPLOT Y= 1.63382*EDF(X,0, 2.43); XMIN=0 XMAX=10 YMIN=0,
    YMAX=0.6 COLOR=1 HEIGHT=2 WIDTH=2
    END
    BEGIN
    FPLOT Y= GDF(X,2.43,1); XMIN=0 XMAX=10 YMIN=0 YMAX=0.4,
    HEIGHT=2 WIDTH=2
    FPLOT Y= 1.1102*GDF(X,2,1.2150); XMIN=0 XMAX=10 YMIN=0,
    YMAX=0.4 COLOR=1 HEIGHT=2 WIDTH=2
    END

    出力は次のようになります。

    図 i 図 ii 図 iii

    一様密度関数 (図 i) をプロポーザルとすると、生成された点のほとんどが受け入れ範囲の外になります。 図 ii (指数) および図 iii (ガンマ) の場合は、ターゲット関数とプロポーザル関数の両方の平均値が同じですが、ガンマ密度関数をプロポーザル関数とすると、定数とプロポーザル関数の積がターゲット関数に近づきます。したがって、プロポーザル関数から生成された点は高い確率でターゲット関数からの標本として受け入れられ、シミュレーション値は理論値 (平均値 = 2.43、分散 = 2.43、 E(X2) =8.3349) に急速に収束します。