Pythonのscipy.optimize.least_squares()で複素数のデータを非線形の複素関数に最小二乗法でフィッティングする。
複素数のデータ、例えばSパラメータを入力にしたとき、あるパラメータを持つ出力関数(これも複素数、かつ普通は非線形)にフィッティングしたいことは良くある。が、scipyの最小二乗法でどうやるのかあまり書いてあるものをみたことがない。
とりあえずあまりスマートじゃないができるようになったので記録として残す。
ポイントは、複素数はoptimize.least_squaresで最小化できないので、いったん実数のラッパーを作ること。
こんな感じ。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
# z2 = a (z1 - b)^2 + c という複素関数とする。
def fit_func(a, b, c, z1):
return a * (z1 - b)**2 + c
def res_func(a, b, c, z1, z2):
return fit_func(a, b, c, z1) - z2
#複素数を使うには実数のラッパーが必要
def res_func_wrap(param, x1, x2):
fwrap = res_func(param[0] + 1j*param[1], param[2] + 1j*param[3], param[4] + 1j*param[5],
x1[:,0]+1j*x1[:,1], x2[:,0]+1j*x2[:,1])
return fwrap.real ** 2 + fwrap.imag ** 2
#推定するa, b, c
a = 1 + 2j
b = -1 + 3j
c = -2 - 1j
#データはn点
n=1000
#z1 と z2を計算(乱数でばらつかせる)
z1 = np.zeros(n, dtype = complex)
z2 = np.zeros(n, dtype = complex)
for i in range(n):
z1[i] = (2 * np.random.rand() - 1) + 1j*(2 * np.random.rand() - 1)
z2[i] = fit_func(a, b, c, z1[i]) + 0.1*( (2 * np.random.rand() - 1) + 1j*(2 * np.random.rand() - 1) )
#実数のアレイに直す。
x1 = np.zeros((n,2))
x2 = np.zeros((n,2))
for i in range(n):
x1[i, 0] = z1[i].real
x1[i, 1] = z1[i].imag
x2[i, 0] = z2[i].real
x2[i, 1] = z2[i].imag
#初期値
initial_values = [0, 0, 0, 0, 0, 0]
res = least_squares(res_func_wrap, initial_values, args = (x1, x2))
print(res.x)
実行すると
[ 1.00113061 2.00065746 -1.00031483 2.9985901 -1.99946214 -1.00700744]
となった。実際は
[1,2,-1,3,-2,-1]
なのでちゃんとパラメータ推定できていそう。ただなんかもっさりしているので、もっといいやり方を誰か知っていたら教えてほしいなあ。
« NKODICE(んこダイス)というのを見て、キヨシ関数みたいだな、、、と思ったのでカシオの高精度計算サイトkeisan.casio.jpに”NKODICE(んこダイス)もどき”を作ってみた。5つがある配列になると??? Scratchでもやってみた。 | トップページ | 待ちに待ったアンデッドガール・マーダーファルス3を読んだ。期待通りめちゃくちゃ面白かった!人外・怪異・物の怪総出演のキャラものでもあり、超燃えるバトルものでもあり、何より青崎有吾さんらしい論理ミステリでトリック(と動機)に驚く! »
「日記・コラム・つぶやき」カテゴリの記事
- 新型コロナウイルス、日本の陽性者数&ワクチン接種者数総計をプロット&中国、韓国、アメリカ、ドイツ、フランス、イギリスの陽性者数もプロット(1/29更新)日本の増加いうほど止まってない…中国はもう5万人で打ち止めということにしたいようだ。9億人感染という話なのに…(2023.01.29)
- JavaScriptの数値計算ライブラリmathjsを使う(11)バーニングシップフラクタルを描いてみる。このココログでも計算できるようにした。(2023.01.31)
- JavaScriptの数値計算ライブラリmathjsを使う(10) リーマンゼータ関数(Riemann Zeta function)を計算、3次元化してPlotlyでぐりぐり動かす。(2023.01.27)
- JavaScriptの数値計算ライブラリmathjsを使う(9) 仏陀のお姿のフラクタル Buddhabrot(ブッダブロ)を描く。このココログでも計算できるようにしてみた。(2023.01.26)
- JavaScriptの数値計算ライブラリmathjsを使う(8) 4段4次のルンゲクッタ法でローレンツ方程式を計算。このココログでもPlotlyで3次元でぐりぐり動かせるようにしてみた。 (2023.01.25)
« NKODICE(んこダイス)というのを見て、キヨシ関数みたいだな、、、と思ったのでカシオの高精度計算サイトkeisan.casio.jpに”NKODICE(んこダイス)もどき”を作ってみた。5つがある配列になると??? Scratchでもやってみた。 | トップページ | 待ちに待ったアンデッドガール・マーダーファルス3を読んだ。期待通りめちゃくちゃ面白かった!人外・怪異・物の怪総出演のキャラものでもあり、超燃えるバトルものでもあり、何より青崎有吾さんらしい論理ミステリでトリック(と動機)に驚く! »
コメント