概要
Pythonで浮動小数点誤差に行き当たった。
コンソール上にはErrorが表示され計算結果はnanになるが、配列を使っていたので、問題がどこで生じているかわかるまでにちょっと手間取った。
解析上の考え
中心が(cx, cy)、半径がrの円を考え、x座標を与えてそれに対する円上のy座標を計算しようとしていたとき。
このような円の方程式は、陰関数では以下のようになる。
これをy ≥ 0で解くと、
問題の発生
中心の座標が(0.5, 0.5)、半径が0.3となるような円が欲しかったので、xの定義域を0.3～0.8として次のようなコードをPythonで書いた。実際はy < 0の部分も必要で、円周上だけでなく内部の座標をランダムに計算し、それらの相関係数を計算するというものだったが、本質的には変わらない。
|
1 2 3 4 5 6 7 8 9 10 |
import numpy as np r = 0.3 cx, cy = 0.5, 0.5 x = np.linspace(0.3, 0.8, 11) y = np.array([cy + np.sqrt(r**2 - (x - cx)**2)]) cor = np.sum((x - x.mean()) * (y - y.mean()) / x.std() / y.std()) print(cor) |
その結果、cor=nanとなってしまう。実行時にsqrtの計算で警告が出ている。
|
1 2 3 |
C:\Users\tomo\Google �h���C�u\IT_and_Mobile\dev\python\trouble\float.py:7: RuntimeWarning: invalid value encountered in sqrt y = np.array([cy + np.sqrt(r**2 - (x - cx)**2)]) nan |
そこで配列x、yを表示させてみると、yの最後の要素がnanになっている。
|
1 2 |
print(x) print(y) |
|
1 2 3 |
[0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 ] [[0.7236068 0.75980762 0.78284271 0.79580399 0.8 0.79580399 0.78284271 0.75980762 0.7236068 0.66583124 nan]] |
さらに、sqrtの中を表示させてみる。
|
1 |
for xx in x: print(r**2, (xx - cx)**2) |
|
1 2 3 4 5 6 7 8 9 10 11 |
0.09 0.04000000000000001 0.09 0.022500000000000006 0.09 0.009999999999999995 0.09 0.0024999999999999988 0.09 0.0 0.09 0.0025000000000000044 0.09 0.010000000000000018 0.09 0.022500000000000006 0.09 0.03999999999999998 0.09 0.0625 0.09 0.09000000000000002 |
最後のところで差し引かれる方の値がごくわずかに大きくなっていて、その結果根号の中が負の値となってnanが発生したらしい。
解決策
linspaceの範囲を円周上の上限・下限から少し内側にすることも考えられるが、今回はnumpy.around()を使った。
numpy.around(a, decimals=0)
これは値aを小数点以下decimalsの桁数で丸めてくれるもので、今回は小数点以下10桁目で丸めてみた。
|
1 2 3 4 5 6 7 8 9 10 |
import numpy as np r = 0.3 cx, cy = 0.5, 0.5 x = np.linspace(0.3, 0.8, 11) y = np.array([cy + np.sqrt(r**2 - np.around(x - cx, 10)**2)]) cor = np.sum((x - x.mean()) * (y - y.mean()) / x.std() / y.std()) print(cor) |
この結果、正常に計算。
|
1 |
-6.625370822144879 |
注意点
- 解析上値が等しくなる（差し引きゼロになる）ような場合でも、浮動小数点ではごくわずかな誤差のため、エラーとなることがある
- numpy.linspace()を使うと上限・下限の間を等分してくれるのでありがたいが、上限・下限自体が他の値と厳密に等しくなければならない場合には、浮動小数点の誤差が生じるため注意が必要