ガウス関数の最小二乗（最小自乗）フィッティング と ガウス窓ミーンシフト

最初に一言

ガウス関数(exp(-x^2))の最小二乗フィッティングの際に、「対数取ってから2次関数近似すればよい」とか嘘だから。

ノイズというのも貴重な情報です。ノイズの分布を考慮せずに空間を安易に変換してはいけません。最小二乗法ではノイズの分布が時刻によって変動しない事を前提としています。

例えば有る時間で信号が1でノイズ幅が1だったとします。べつの時刻で信号が100でノイズ幅が1だとします。これは最小二乗法でフィッティングしてＯＫです。

ここで対数を取るとざっくりlog(2)-log(1)=0.69とlog(101)-log(100)=0.0099となって小さい信号の上にのったノイズにひっぱられて、大きな信号の上のノイズを見逃します。

同様に、多くの場合で最尤推定!=最小二乗です。

あーーー、すっきりした。　

では本題

ということでガウス関数の最小二乗をやろうと思ってググってたんですが、なんか良いサンプルが見つからないんですね。で、頑張って自分で微分して、解を求めようとしたんですが、どーも計算できない。（微分・積分苦手）

最尤推定のときは解析解が求まるけど、最小二乗だと解析解がもとまらないのかな？という事で普通はニュートン法に逃げる訳だけれども、微分式を見るとなんか既視感。

g(t)をガウスモデル関数とする（なお、Σ時間tでの積算）

振幅aで微分: Σa*g(t)*g(t) =Σg(t)*i(t)

平均uで微分: Σa*g(t)*g(t)*t =Σg(t)*i(t)*t

標準偏差sで微分:Σa*g(t)*g(t)*t*t =Σg(t)*i(t)*t*t

モーメントの形に成っちゃってるけれども、ざっくりと言うと。

入力信号にモデル関数を掛けた平均値が、

モデル関数にモデル関数をかけたものの平均値になるようにパラメータを決めなさい。

ってこったね。

ん？ここではガウス窓で見てるけれども、矩形窓だと考えて平均値のみを考えるとミーンシフトそのものだよね？

1.適当な矩形範囲で重心を求める

2.上で求めた重心を中心にした矩形範囲で重心を求め直す

3.落ち着くまで2を繰り返す

ということで、ミーンシフトとして解釈して実装してみた。ただし、全体を最小二乗マッチングする必要あるので窓範囲=入力信号範囲。そのため、窓関数の重心が中心にくるとは限らない状況。

なんでモデルの重心も意識しながらコーディングする必要があるのがオリジナルのミーンシフトと違うという解釈。

結局は原始的な山登り法なんじゃね？という疑惑がふつふつと自分の中でわき上がるけれども、今後非対称重みのミーンシフトを実装する自信がついたのでよしとする。

以下rubyのサンプルコード（標準偏差もうごいてるけど真面目に考えてない)

#!/usr/bin/ruby
# -*- coding: utf-8 -*-
include Math

#テスト用入力信号
input=Array.new
100.times do |i|
signal = exp(-((i.to_f-75)/30.0)**2/2.0)*4.0;
signal += 0.1 * rand #-0.5;
input.push(signal)
end

#ガウス関数
##t=time u= center s=stddev
def gauss(t,u,s)
return exp(-((t-u)/s)**2/2)
end

#初期値を適当に設定(対数とって最小二乗でもいいよ)
m0=0.0
m1=0.0
m2=0.0
input.size.times do |t|
m0+=input[t]
m1+=input[t] * t
m2+=input[t] * t**2
end
u= m1/m0;
s= sqrt(m2/m0 - u**2);

#振幅を求める
def detectA(input,u,s)
gg=0.0
gi=0.0
input.size.times do |t|
g=gauss(t,u,s);
gg += g**2;
gi += g*input[t]
end
return gi/gg;
end

#平均と標準偏差更新（ついでに振幅)
def detectUS(input,a,u,s)
mg0=mg1=mg2=mi0=mi1=mi2=0.0
input.size.times do |t|
g=a*gauss(t,u,s);
gi=g * input[t]
gg= g* g
mi0+=gi
mi1+=gi * t
mi2+=gi * t**2

mg0+=gg
mg1+=gg * t
mg2+=gg * t**2


end
ui=mi1/mi0 #実際の重心候補
ug=mg1/mg0 #理論上の重心 (データ打ち切りの影響で中心!=重心）
vi=mi2/mi0 - ui**2 #実際の分散候補
vg=mg2/mg0 - ug**2 #理論上の分散
a=detectA(input,u,s)
return [a,u+ui-ug,sqrt(s**2 + vi-vg) ]
end
a=detectA(input,u,s)

#更新処理
100.times do |i|
d=detectUS(input,a,u,s)
a=d[0]
u=d[1]
s=d[2]
puts "#{a}\t#{u}\t#{s}";
end

#結果出力
open("o","w")do |f|
input.size.times do |t|
f.puts "#{input[t]}\t#{ a * gauss(t,u,s)}"
end
end