2016-06-15 8 views
1

データ生成関数はかなり面倒なので、私はできるだけ明確にしようとします。もしそうでないとコメントしてください。Scipyのクワッドが「きれいな」機能で失敗する

私はいくつかの変数に関数を持ち、もう一方を一定に保ちながら積分しようとしています。しかし、そのセカンダリ変数の値が異なると、統合プロセスは全く異なる結果になります。

ここに私のサンプルコード(。unfort再現可能ではない)です。

fig, ax = plt.subplots() 
tPrimeGrid = [0.16, 0.18, 0.22] 
from scipy import integrate 
for ttPrime in tPrimeGrid: 
    # compute grid 
    lowerBar = continuousTime.getPLowerBarAnalytical(ttPrime, Param.S, Param) 
    upperBar = Param.r 

    # integrate function 
    h = integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, ttPrime, Param), full_output=True) 
    print('tPrime {} value {} erorr {}'.format(ttPrime, h[0], h[1])) 
    # plot function to be evaluated 
    pGrid = np.linspace(lowerBar, upperBar, 100) 
    plt.plot(pGrid, [myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid], label='t {} (h: {})'.format(ttPrime, h[0])) 
plt.legend() 

quadの出力

tPrime 0.16 value 6.63310536371 erorr 0.000345564616621 
tPrime 0.18 value 5.93645658492 erorr 0.00045956820487 
tPrime 0.22 value 0.359208009237 erorr 3.98801002485e-15 

エラーがうまく寛容の私のレベル内のすべてのですが、値は多くのことをジャンプさ。ここではグラフです:

plot of functions

したがって、これらの機能は非常によく似た形状を有しています。肉眼で見ると、積分の違いは実際には意味をなさない。そこで、私はそれらを「手動で」計算しました。

(np.array([myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid])).sum()*(pGrid[1]-pGrid[0]) 
0.35538925924926557 

小さな値が実際には正しいことを意味します。だからここ

integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True) 
Out[19]: 
(6.634157704675579, 
0.004721148834250418, 
{'alist': array([ 1.99994895, 1.78826738, 1.86060326, 1.79090489, 1.93030163, 
      1.72120652, 1.96515082, 1.75605571, 1.98257541, 1.7734803 , 
      1.9912877 , 1.7821926 , 1.99564385, 1.78872682, 1.99782193, 
      1.78654874, 1.99940018, 1.78763778, 1.99945548, 1.99891096, 
      1.78845456, 1.99972774, 1.99918322, 1.78831843, 1.99988939, 
      1.99931935, 1.7881823 , 1.99999149, 1.99942145, 1.7882844 , 
      1.9998979 , 1.9999447 , 1.99940443, 1.78825036, 1.99996597, 
      1.99986387, 1.99991492, 1.99995746, 1.99938742, 1.78827589, 
      1.99993194, 1.99990641, 1.99998298, 1.99988089, 1.99995321, 
      1.99939592, 1.78827164, 1.99994044, 1.99995108, 1.99940231]), 
    'blist': array([ 1.99995108, 1.78827164, 1.93030163, 1.86060326, 1.96515082, 
      1.75605571, 1.98257541, 1.7734803 , 1.9912877 , 1.7821926 , 
      1.99564385, 1.78654874, 1.99782193, 1.79090489, 1.99891096, 
      1.78763778, 1.99940231, 1.7881823 , 1.99972774, 1.99918322, 
      1.78872682, 1.99986387, 1.99931935, 1.78845456, 1.9998979 , 
      1.99938742, 1.78825036, 2.  , 1.99945548, 1.78831843, 
      1.99990641, 1.99994895, 1.99942145, 1.78826738, 1.99998298, 
      1.99988089, 1.99993194, 1.99996597, 1.99939592, 1.7882844 , 
      1.99994044, 1.99991492, 1.99999149, 1.99988939, 1.99995746, 
      1.99940018, 1.78827589, 1.9999447 , 1.99995321, 1.99940443]), 
    'elist': array([ 4.61134496e-03, 4.14135579e-06, 1.27804393e-15, 
      1.55808958e-15, 5.54429458e-16, 0.00000000e+00, 
      2.90617202e-16, 0.00000000e+00, 2.27720143e-15, 
      0.00000000e+00, 3.91514838e-15, 0.00000000e+00, 
      2.15987477e-14, 5.40749179e-17, 1.19682144e-12, 
      0.00000000e+00, 1.03521880e-04, 0.00000000e+00, 
      2.48275100e-08, 6.85735811e-13, 6.78454062e-18, 
      8.65204618e-09, 1.64268148e-13, 3.39437556e-18, 
      4.65487056e-08, 1.12644353e-13, 0.00000000e+00, 
      6.19443638e-08, 4.08066769e-09, 8.48813317e-19, 
      5.99147851e-07, 1.21478039e-06, 9.39741081e-10, 
      0.00000000e+00, 6.32087778e-14, 2.19912539e-10, 
      6.97448944e-09, 1.85114515e-13, 1.28558440e-13, 
      2.12217047e-19, 8.89451362e-08, 8.45167083e-09, 
      1.25621961e-14, 2.59393721e-08, 2.45738052e-15, 
      1.19106919e-12, 1.06110581e-19, 4.91812027e-08, 
      2.72831861e-15, 3.06819516e-12]), 
    'iord': array([ 1, 17, 50, 49, 48, 47, 46, 45, 44, 43, 42, 29, 40, 39, 38, 23, 35, 
     11, 34, 3, 7, 21, 30, 18, 12, 8, 6, 0, 0, 0, 0, 0, 0, 0, 
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int32), 
    'last': 50, 
    'neval': 2079, 
    'rlist': array([ 1.04205921e-01, 6.52426361e-06, 1.15115963e-01, 
      1.40340233e-01, 4.99385660e-02, 0.00000000e+00, 
      2.33230318e-02, 0.00000000e+00, 1.12779305e-02, 
      0.00000000e+00, 5.54633015e-03, 0.00000000e+00, 
      2.75040018e-03, 4.87063561e-03, 1.36955727e-03, 
      0.00000000e+00, 8.85474806e-03, 0.00000000e+00, 
      1.73731989e+00, 3.41803845e-04, 6.11097092e-04, 
      1.73678061e+00, 1.70814248e-04, 3.05738170e-04, 
      2.00530044e-01, 8.53852175e-05, 0.00000000e+00, 
      5.29681716e-06, 1.51991445e-01, 7.64543067e-05, 
      2.17985628e-01, 2.00512965e-01, 7.26773425e-02, 
      0.00000000e+00, 1.06564581e-05, 3.34545761e-01, 
      5.59012296e-01, 5.32838374e-06, 1.06721256e-05, 
      1.91148122e-05, 3.34512038e-01, 2.38773007e-01, 
      5.32807434e-06, 1.85664300e-01, 2.66423055e-06, 
      5.33597722e-06, 9.55759147e-06, 1.85647516e-01, 
      1.33212494e-06, 8.93844378e-03])}, 
'The maximum number of subdivisions (50) has been achieved.\n If increasing the limit yields no improvement it is advised to analyze \n the integrand in order to determine the difficulties. If the position of a \n local difficulty can be determined (singularity, discontinuity) one will \n probably gain from splitting up the interval and calling the integrator \n on the subranges. Perhaps a special-purpose integrator should be used.') 

を私の関連質問です:したがって、私も悪いものの一つはここをquad()の完全な出力を追加します

  1. はどのようにそれができるような同様の形状で、いくつかのパラメータは動作し、いくつかは動作しません。これらの失敗を将来どのように「予測」することができますか?
  2. 完全な出力(以下)は、分割の最大数が達成されたことを示しています。エラーが潜在的に見積もられていないとは言えません。それは意味する必要がありますか? 5.9-0.00040.3の近くにありません。
  3. 制限値が助けにならなかったので(下記参照)、潜在的な選択肢は何ですか?この機能をどのように統合できますか?

出力(良くない)私はlimitを増やしてみましたが、それは私だけ次を与えた:

integrate.quad(expectedUtility, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True, limit=500) 
Out[24]: 
(6.633112814769514, 
4.74743687572826e-06, 
[...] 
'The occurrence of roundoff error is detected, which prevents \n the requested tolerance from being achieved. The error may be \n underestimated.') 

答えて

2

quadfull_outputがTrueで、戻り値のため、多少ファンキーな規則があります。 quadが問題を検出しない場合は、y, abserr, infodictを返します。 quadが何らかの形式の障害を検出すると、y, abserr, infodict, messageを返します。 infodictには、失敗を示す単純なstatusフィールドが含まれていないため、4番目の戻り値の有無を確認する必要があります。存在する場合は、問題を説明する文字列です。 (。あなたのコードでは、あなたがif len(h) > 3: ...handle the failure...のようなものが追加される場合があります)悪い例あなたの完全な出力では、あなたはmessageこれを含む文字列であることがわかります。

The maximum number of subdivisions (50) has been achieved. 
If increasing the limit yields no improvement it is advised to analyze 
the integrand in order to determine the difficulties. If the position of a 
local difficulty can be determined (singularity, discontinuity) one will 
probably gain from splitting up the interval and calling the integrator 
on the subranges. Perhaps a special-purpose integrator should be used. 

これは数値積分が失敗したことを意味します。 integrand myFuncについて詳しく知ることなく、それがなぜ失敗したのかは分かりません。

関連する問題