Pythonで、複数の座標情報から、その座標を通る三角関数や指数・対数関数の式あるいは近似式を導出したい

Question

数学寄りの質問で大変申し訳ないのですが、
Pythonで、複数の座標情報（座標情報は、いくつでも取得できることを想定。もちろん、利用する座標情報はなるべく少ない方が良い）を元に、その座標を通る三角関数や指数・対数関数の式あるいは近似式を導出したいと考えております（その式の係数を求めることになります）。なお、三角関数、指数・対数関数のどの関数の式であるかがわかっている前提とします。

今のところ、うまく導き出す方法を考えられておらず、
・座標を用いて計算を行い、導き出す
・Pythonで利用出来る関数に、導出してくれる関数があり、それを用いることで導き出す
・wolframのようなサイトがあって、そこのAPIに導出してくれる関数があり、それを用いることで導き出す
・例えば正弦波なら、y = asin(bx) + cの式の、a, b, cのとりうる範囲を限定し、二分探索などで、座標が一致しているかチェックすることで導き出す
という４つのぼんやりした方法しか思いつきません。

ネット上で関連した記事を探してみたのですが、Excelのソルバーや物理学（正弦波、余弦波、正接波など）などの記事しか見つかりませんでした。

なにか他に良い、式の導出（式の係数を求める）方法はありませんでしょうか？
プログラム関連のQ&Aサイトにこういった質問を投稿するのはお門違いな気もするのですが、どなたかご教授頂けないでしょうか？

Comments

[swordone]

与えられた座標を通る(または近似する)三角関数なり指数・対数関数なりの係数を求めるという意味に見えるのですが，違うのですか？

[tukejonnyBot]

はい。そうなります。
説明不足でした。申し訳ありません・・・

Answer 1

最小二乗法が使いやすいと思います．
いくつか追加のモジュールが必要です．
Macで試したのですが，NumPyとMatplotLibとSciPyをインストールすれば良いと思います．

$ pip3 install numpy
$ pip3 install matplotlib
$ pip3 install scipy

その上で，以下のコードを実行してみてください．
テストに使ったデータは y = x ** 2 + 1 の座標です．
２次関数の一般形 y = a * x ** 2 + b * x + c を定義して，
そのパラメータを推定しています．
a_fit，b_fit，c_fitにそれぞれ推定したパラメータが格納されます．
三角関数，指数・対数関数などでも利用できると思うので，試してみてください．

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

x = np.array([0, 1, 2, 3, 4, 5, 6, 7])
y = np.array([1, 2, 5, 10, 17, 26, 37, 50])

def fit_func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    c = parameter[2]
    residual = y-(a*x**2+b*x+c)
    return residual

parameter0 = [0.,0.,0.]
result = optimize.leastsq(fit_func,parameter0,args=(x,y))
a_fit=result[0][0]
b_fit=result[0][1]
c_fit=result[0][2]
print(a_fit,b_fit,c_fit)

#PLot
plt.figure(figsize=(8,5))
plt.plot(x,y,'bo', label='Exp.')
plt.plot(x,a_fit*x**2+b_fit*x+c_fit,'k-', label='fitted parabora', linewidth=10, alpha=0.3)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.legend(loc='best',fancybox=True, shadow=True)
plt.grid(True)
plt.show()

参考
フィッティング３(最小二乗法，optimizeの利用) — Mechanical Design Lab. of TUMSAT 1.0 documentation

Comments

[tukejonnyBot]

回答ありがとうございます！！
最小二乗法という方法があるのですね・・・

勉強不足で大変申し訳ないのですが、一つわからないところがあります。
optimize.leastsq()の返り値が二次元配列になっているのですが、値はどのように格納されているのでしょうか？
単純にパラメータ１つ目からresult[0][1], result[0][2], ..., result[1][0], result[1][1],...となっているのでしょうか？

実際に三角関数(sin)のグラフを描画してみました。描画されているグラフはあっていそうなのですが、パラメータが少しおかしいように思えます

```Python
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from math import*

""" test data """
x = np.array([r for r in range(1, 1000)])
y = np.array([sin(pi * (r) / 180.0) for r in range(1, 1000)])
""" test data """

def calcSin(angle):
rad = pi * angle / 180.0
for r in range(0, len(angle)):
angle[r] = sin(rad[r])
return(angle)

def fit_func(parameter, x, y):
a = parameter[0]
b = parameter[1]
c = parameter[2]

residual = y-(a*calcSin(b * x)+c)
return(residual)

parameter0 = [0.,0.,0.]
result = optimize.leastsq(fit_func, parameter0, args=(x, y))
print "result is " + str(result)
a_fit = result[0][0]
b_fit = result[0][1]
c_fit = result[0][2]

print "y = (" + str(a_fit) + ")sin((" + str(b_fit) + ")x) + " + str(c_fit)

#PLot
plt.figure(figsize=(8,5))
plt.plot(x,y,'bo', label='Exp.')
plt.plot(x,a_fit*calcSin(b_fit * x)+c_fit,'k-', label='fitted parabora', linewidth=10, alpha=0.3)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.legend(loc='best',fancybox=True, shadow=True)
plt.grid(True)
plt.show()
```

出力結果（グラフは除く）
y = (0.0)sin((0.0)x) + 0.0478855562909

これだとxの値に関わらず、yが0.0478855562909となるように見えてしまいます。
出力結果として、これは正しいのでしょうか？

[KenTerada]

scipy.optimize.leastsqの返り値がresultだったとき，
これは２次元配列になっていて，
result[0]が推定パラメータ，result[1]がパラメータの誤差だそうです．
result[0]ばかり気にしていたので，result[1]が真に何を意味するのかはわかりません．

Answer 2

コメントだとコードが見にくくなるので，別の所へ．
確かに，色々やってみましたが上手く行きませんでした．
入力データを変えた所上手く動いたのですが，なぜだかわからずじまいです．
もしかしたらどこかで，浮動小数の情報落ち，桁落ちなどで
値の更新がされなくなってしまったのかもしれません．
（最小二乗法は，計算を何度も繰り返して最小になる値を求めるので）

pylabのsin関数は配列に対しても行えると聞いたので，使っています．
よくよく考えたら，yの式をmyfitか何かで宣言しておいて，
mydiffか何かで入力データとの差を取れば良いかなと思いました．

また，初期値も[0.,0.,0.]でやるとエラーになることを確認しています．
最小二乗法では最小値を推定する計算があるのですが，
その方法を使う時にいくつか制限があるということです．
この場合は１番目の値が0だとだめなようです．
（ヤコビ行列が陽，と私にもさっぱりの定義でした）

import matplotlib.pyplot as plt
from scipy import optimize
from math import *
import pylab

""" test data """
x = [(pi * (r) / 180.0) for r in range(1, 1000)]
y = pylab.sin(x)
""" test data """

def myfit(p, x):
	return(p[0]*pylab.sin(p[1]*x)+p[2])

def mydiff(p, x, y):
	return(y-myfit(p,x))

p0 = [1.0,1.0,1.0]
result = optimize.leastsq(mydiff, p0, args=(x,y))
pr = result[0]
print pr

#Plot
plt.figure(figsize=(8,5))
plt.plot(x,y,'bo', label='Exp.')
plt.plot(x,myfit(pr,x),'k-', label='fitted curb', linewidth=10, alpha=0.3)
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.legend(loc='best',fancybox=True, shadow=True)
plt.grid(True)
plt.show()