正確な警告はscipy.specialで

Question

私はscipy.special.hyp2f1()を使用して各ステップで超幾何関数の計算を必要とするMCMCサンプラーを実行しています。私のグリッド上の特定のポイントで

（私は気にしない）超幾何関数の解決策は非常に不安定であり、scipyのダウンロードは警告を出力します。

Warning! You should check the accuracy 

これはかなり面倒ですし、サンプルの1000を超えます私のルーチンを良くするかもしれない。

私はspecial.errprint(0)を使用してみましたが、warningsモジュールと-W ignoreフラグの両方を使用してPythonのすべての警告を無効にしました。

（別のファイルから呼び出された）問題の機能は

from numpy import pi, hypot, real, imag 
import scipy.special as special 

def deflection_angle(p, (x1, x2)): 

    # Find the normalisation constant 
    norm = (p.f * p.m * (p.r0 ** (t - 2.0))/pi) ** (1.0/t) 

    # Define the complex plane 
    z = x1 + 1j * x2 

    # Define the radial coordinates 
    r = hypot(x1, x2) 

    # Truncate the radial coordinates 
    r_ = r * (r < p.r0).astype('float') + p.r0 * (r >= p.r0).astype('float') 

    # Calculate the radial part 
    radial = (norm ** 2/(p.f * z)) * ((norm/r_) ** (t - 2)) 

    # Calculate the angular part 
    h1, h2, h3 = 0.5, 1.0 - t/2.0, 2.0 - t/2.0 
    h4 = ((1 - p.f ** 2)/p.f ** 2) * (r_/z) ** 2 
    special.errprint(0) 
    angular = special.hyp2f1(h1, h2, h3, h4) 

    # Assemble the deflection angle 
    alpha = (- radial * angular).conjugate() 

    # Separate real and imaginary parts 
    return real(alpha), imag(alpha)` 

Answer 1

下回る残念ながら、hyp2f1はパラメータ空間のいくつかの非自明な領域を超える計算するために悪名高く難しいです。多くの実装では、間違った結果や間違った結果が生じることがあります。 Scipy.specialは少なくともコンバージェンスを監視するのは難しいです。代わりに、任意の精度の実装をusrにすることもできます。 mpmath。しかし、これらは確かにかなり遅いので、MCMCユーザーは注意してください。

編集：これはscipyのバージョンに依存します。 @ wrwrwrの例をscipy 0.13.3で試してみました。見た目を再現しています。errprintのステータスに関係なく、「警告！正確さをチェックしてください」と表示されます。しかし、dev版でも同じことをやっています。

In [12]: errprint(True) 
Out[12]: 0 

In [13]: hyp2f1(0.5, 2/3., 1.5, 0.09j+0.75j) 
/home/br/virtualenvs/scipy_py27/bin/ipython:1: SpecialFunctionWarning: scipy.special/chyp2f1: loss of precision 
    #!/home/br/virtualenvs/scipy_py27/bin/python 
Out[13]: (0.93934867949609357+0.15593972567482395j) 

In [14]: errprint(False) 
Out[14]: 1 

In [15]: hyp2f1(0.5, 2/3., 1.5, 0.09j+0.75j) 
Out[15]: (0.93934867949609357+0.15593972567482395j) 

だから、明らかに2013年と今の間にある時点で修正されています。あなたはscipyのバージョンをアップグレードしたいかもしれません。