問題タブ [numerical-stability]

私は機能を持っていました：

見たことがない人にとってはx0 + (x1-x0) * alpha、後者はそれを保証しないため、これよりも望ましいlerp(1.0f, x0, x1) == x1.

ここで、lerp関数に追加のプロパティが必要です: I'd like lerp(alpha, x0, x1) == lerp(1-alpha, x1, x0). （理由については、これはより複雑な関数のおもちゃの例です。）私が思いついた、うまくいくように見える解決策は

この二重減算には、0 付近と 1 付近を丸める効果がalpha = std::nextafter(0)あり1 - alpha == 1ます1 - (1-alpha) == 0。私が知る限り、それは常に真実です1.0f - x == 1.0f - (1.0f - (1.0f - x))。という効果もあるそうですw0 + w1 == 1.0f。

質問:

1. これは合理的なアプローチですか？
2. 自分のコンパイラが自分のやりたいことをやってくれると信頼できますか? 特に、Windows では、部分的な結果に対してより高い精度を使用することがあることを知っています。また、コンパイラが代数を実行できることも知っています。代数的に明らかに 1-(1-x)==x です。

これは、Clang、VisualStudio、および gcc を使用した C++11 です。

線形の生死過程の遷移確率を数値的に評価したい

どこで 二項係数と

ほとんどのパラメーターの組み合わせについて、(対数と Kahan-Neumaier 加算アルゴリズムを使用して) 許容可能な数値誤差で評価できます。

加数の符号が交互になり、数値誤差が合計を支配する場合に問題が発生します (この場合、条件数は無限大になる傾向があります)。これは次の場合に発生します。

たとえば、評価に問題がありますp(1000, 2158, 72.78045, 0.02, 0.01)。0 のはずですが、非常に大きな値log(p) ≈ 99.05811が得られました。これは、確率としては不可能です。

さまざまな方法で合計をリファクタリングし、Zhu-Hayesなどのさまざまな「正確な」合計アルゴリズムを使用してみました。私はいつもほぼ同じ間違った値を取得しているため、問題は数値を合計する方法ではなく、各加数の内部表現にあると思います。

二項係数のため、値は簡単にオーバーフローします。各 (絶対) 要素を最小正規数と 1 の間の合計に保つために、線形変換を試みました。

私は今行き止まりにあり、先に進む方法がわかりません。任意精度の算術ライブラリを使用できますが、マルコフ連鎖モンテカルロ アプリケーションには計算コストが高すぎます。

IEEE-754 double に十分な精度で部分和を格納できない場合に、そのような和を評価する適切な方法またはトリックはありますか?

これは、Kahan 加算アルゴリズムを使用して最大値と合計で値を再スケーリングするだけの基本的な作業例です。明らかに、ほとんどの値は Float64 ではサブノーマルになります。

Python3 を使用してオイラー法でこの微分方程式を解こうとしています。

Wolfram Alpha によると、これが正しい方程式のプロットです。

繰り返しになりますが、Wolfram Alpha によると、この場合、間隔の終わりまでにわかるように、従来のオイラー法は安定していないはずです。

ただし、私の実装では、オイラー法は安定した結果を提供しますが、これは奇妙なことです。何らかの理由で私の実装が間違っているのだろうか。それにもかかわらず、私はエラーを見つけることができません。

近似値と関数の分析出力を比較するいくつかの点とプロットを生成しました。青は対照群としての分析結果。赤で、私の実装の出力:

それが私のコードです：

助けてくれてありがとう。

================================================== ==========

アップデート

@SourabhBhat の助けを借りて、実装が実際に正しいことを確認できました。それは確かに、不安定さを生み出していました。ステップサイズを大きくするだけでなく、ズームインしてそれが起こっていることを確認する必要がありました.

次の図は、それ自体を物語っています (ステップ サイズ 0.22)。