2

私はモンテカルロアルゴリズムを書いています。このアルゴリズムでは、ある時点で確率変数で除算する必要があります。より正確には、確率変数は差分商のステップ幅として使用されるため、実際には最初に変数を乗算してから、この式の局所的な線形関数から再度除算します。好き

double f(double);

std::tr1::variate_generator<std::tr1::mt19937, std::tr1::normal_distribution<> >
  r( std::tr1::mt19937(time(NULL)),
     std::tr1::normal_distribution<>(0) );

double h = r();
double a = ( f(x+h) - f(x) ) / h;

これはほとんどの場合正常に機能しますが、。の場合は失敗しh=0ます。数学的には、これは問題ではありません。正規分布の確率変数の有限（または実際には可算）選択では、すべてが確率1で非ゼロになるためです。しかし、デジタル実装では、h==0約2³²の関数呼び出しが発生します。 （宇宙よりも周期が長いメルセンヌツイスターに関係なく、通常longのsを出力します！）。

この問題を回避するのは非常に簡単です。

double h = r();
while (h==0) h=r();

しかし、私はこれを特にエレガントだとは思いません。より良い方法はありますか？

私が評価している関数は、実際には単純なℝ->ℝのようなものfではなく、ℝⁿ変数を数値積分しながらℝᵐ変数の勾配を計算するℝᵐxℝⁿ->ℝです。関数全体に、予測できない（しかし「コヒーレントな」）ノイズが重なっており、特定の（しかし未知の）優れた周波数がある場合もあります。これを固定値で試してみると、問題が発生しますh。

4

5 に答える 5

3

あなたのやり方は十分にエレガントに見えますが、少し違うかもしれません：

do {
    h = r();
} while (h == 0.0);
于 2011-06-18T22:19:34.887 に答える
3

2つの正規分布確率変数の比率は、コーシー分布です。コーシー分布は、分散が無限大の厄介な分布の1つです。確かに非常に厄介です。コーシー分布は、モンテカルロ実験を混乱させます。

2つの確率変数の比率が計算される多くの場合、分母は正常ではありません。人々はしばしば正規分布を使用して、この非正規分布の確率変数を近似します。

  • 正規分布は通常、操作が非常に簡単です。
  • 通常、そのような素晴らしい数学的特性を持っています、
  • 通常の仮定は多かれ少なかれ正しいように見えます、そして
  • 実際の分布はクマです。

距離で割っているとします。距離は、定義上、半正定値であり、確率変数として正定値であることがよくあります。したがって、すぐに距離を正規分布させることはできません。それにもかかわらず、平均が標準偏差よりもはるかに大きい場合、人々は距離の正規分布を想定することがよくあります。この通常の仮定が行われるとき、あなたはそれらの非現実的な値から保護する必要があります。簡単な解決策の1つは、切断正規分布です。

于 2011-06-18T22:35:25.597 に答える
0

計算しようとしているものによっては、おそらく次のようなものが機能します。

double h = r();
double a;
if (h != 0)
    a = ( f(x+h) - f(x) ) / h;
else
    a = 0;

fが線形関数である場合、これは（私が思うに？）h=0で連続のままである必要があります。

代わりに、ブランチのコストを回避するために、ゼロ除算の例外をトラップすることを検討することもできます。これはパフォーマンスに悪影響を与える場合とそうでない場合があることに注意してください-両方の方法でベンチマークを行ってください！

-fnon-call-exceptionsLinuxでは、ゼロによるゼロ除算の可能性を含むファイルを作成し、SIGFPEハンドラーをインストールする必要があります。

struct fp_exception { };

void sigfpe(int) {
  signal(SIGFPE, sigfpe);
  throw fp_exception();
}

void setup() {
  signal(SIGFPE, sigfpe);
}

// Later...
    try {
        run_one_monte_carlo_trial();
    } catch (fp_exception &) {
        // skip this trial
    }

Windowsでは、SEHを使用します。

__try 
{ 
    run_one_monte_carlo_trial();
} 
__except(GetExceptionCode() == EXCEPTION_INT_DIVIDE_BY_ZERO ? 
         EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH)
{ 
    // skip this trial
}

これには、高速パスへの影響が少ない可能性があるという利点があります。例外ハンドラレコードの調整があるかもしれませんが、ブランチはありません。Linuxでは、コンパイラがforのより保守的なコードを生成するため、パフォーマンスがわずかに低下する可能性があります-fnon-call-exceptions。下でコンパイルされたコードが-fnon-call-exceptions自動（スタック）C ++オブジェクトを割り当てない場合、これが問題になる可能性は低くなります。これにより、ゼロによる除算が非常に高額になる場合もあります。

于 2011-06-18T22:14:35.523 に答える
0

関数 (f(x+h)-f(x))/h には h->0 の制限があるため、h==0 に遭遇した場合はその制限を使用する必要があります。制限は f'(x) になるので、微分がわかっている場合はそれを使用できます。

実際に行っているのは、正規分布に近似する離散点の数を作成することであり、これが分布にとって十分である場合は、実際に値が 0 にならないように作成します。

于 2011-06-18T22:29:06.870 に答える
0

正規分布を維持したい場合は、0 を除外するか、以前は発生しなかった新しい値に 0 を割り当てる必要があります。コンピューター サイエンスの限られた範囲では 2 番目の方法はおそらく不可能であるため、1 番目の方法が唯一の選択肢です。

于 2011-06-18T22:13:36.613 に答える