非線形最小二乗によるコールコールパラメータの求め方

質問

非線形最小二乗法によりコールコールパラメータを求めたいのですがよくわかりません。
詳細は以下の通りです。

コールコールの式ρ(iω)=ρ(0)*(1-m*(1-(1/(1+(iωτ)C))))の、複素比抵抗ρ(iω)が測定データであり、ωが測定周波数でこれが既知の値となっています。
ここから
ρ(0)：直流比抵抗
m：充電率
τ：時定数
ｃ：周波数依存係数
の４つのパラメータを非線形最小二乗法（修正マルカート法など）を用いてρ(iω)=ρ(0)*(1-m*(1-(1/(1+(iωτ)C))))にフィッティングさせて求めたいのですが上手くいきません。コールコールの式のように複素数の式であり、さらにiにc乗がかかっている場合はどのように考えたらいいのでしょうか？例えば
誤差二乗和（χ2）→実部と虚部に分けて考え誤差二乗和を求めるときに足すのか？
ヤコビアン行列→複雑な式を無理やり各パラメータで偏微分するのか？
などです。（上記に挙げた考え方自体間違っているかもしれません。。。）

どなたか複素インピーダンスなど複素数の式を数値計算する時の考え方など、経験がある方がいらっしゃればアドバイスして頂きたいです。宜しく御願い致します。

回答No.2

vq100mg

1 / ( 1 + (iωτ)^C ) という形を含む類似の緩和曲線フィッティングの経験があります。

私の場合も４パラメータでした。　評価は最小二乗で、最急降下法と準ニュートン法を試しました。　各パラメータを微小変化させ、それぞれに対する残差変化(感度)に基づいて探索方向を反復修正していく、すなわち各パラメータを動かす比率(加重)を逐次更新していき、最小二乗誤差点に向かわせるというものです。　偏微分は概念的なものであり、その値は数値的に得るだけです。　僅かに異なる数値を代入して傾きを求めるのです。　解析的微分の手計算は一切ありません。

複素数の間の距離という一般論からすれば、Σ ( (実部同士の差)^2 + (虚部同士の差)^2 ) を評価関数とするのは妥当でしょう。　しかしもし、あなたの式の m の値がかなり小さいなら、虚数部のみを使用した方が良いかもしれません。　この場合、特徴は虚数部に十分出ており、大きなオフセットを持ちこむ実数部は感度的にむしろ邪魔な存在でしょう。

ところで周波数点は、対数的に選んでいらっしゃるでしょうか。　ご提示の曲線は対数周波数で見て対称的なものです。　言うまでも無いと思いますが、例えばリニアに周波数点を選ぶなど、対数周波数で見て評価点密度が低域に偏ると、適正なフィッティングが出来ないかもしれません。

お礼

vq100mgさん

回答ありがとうございます。各パラメータのオーダーですが、だいたい
ρ(0)=50～500
m=0.1～2.0
τ=0～1.0
c=0～1.0
くらいだと思います。ω（角周波数）は対数的に選んでいます。
評価関数ですが虚数部のみを使用して試してみます。フィッティングさせるコールコールの式を変形させたりしないとダメと思っていましたが
数値計算上のテクニックで求めることができるのですね。回答を参考にさせて頂き数値計算の勉強をもう一度してみます。

回答No.1

usokoku

やりましたけど、具体的テクニックを書くと個人が特定できる(あんな馬鹿げた内容をするのは私ぐらい)ので、一般論だけ

＞誤差二乗和
計算値が発散しやすいので、残差絶対値和を使用
＞偏微分するのか？
偏微分で求めた値を使うと大体発散してくれるので使わない。

初期値のあたりを、たとえば、ピポットをひとつ選んで、100ぐらいに分割した範囲で残差を計算して、もっとも減った値を選定する
という、ORに近い考え方で計算して行くと、発散しにくいです。
どの値をピポットとするかは、上下１０個くらいを計算して、作図して、形の崩れない点を選びます。収束が早い点では、グラフの形が崩れて、カイノス問題とか、鞍部問題とかで、地獄を見ます。

お礼

usokokuさん

回答ありがとうございます。一度マルカート法でプログラムを作成したことはあるのですが、実は数値計算はあんまりやった事ないので完璧に理解していません。回答を参考にさせて頂き、数値計算もう一度勉強して見たいと思います。