最小二乗法の重みづけについて教えてください。

質問

簡単のため、説明変数tと、目的変数xが、共に実数（スカラー）とします。
また、フィッティング関数
F=F(t,a,b,c)
も、簡単のため3変数または4変数のスカラー値関数とし、フィッティングパラメータa,b,cも実数(スカラー)とします。また、Fがフィッティングパラメータを2つしか持たない場合（Fが3変数の場合）には、
F(t,a,b)と読み替えて考えることにします。

また、データ、即ち説明変数と目的変数の実測値の組　(t_{i},y_{i})がn個あるとする。

また、以下の４種類の評価関数を考えます。

A:所謂2ノルム
A(a,b,c)=Σ|{y}_{i}-F({t}_{i},a,b,c)|^2

B:n個の正数w_{i}を用いて、重みづけ
B(a,b,c)=Σ(w_{i}|{y}_{i}-F({t}_{i},a,b,c)|^2)

C:単調（非退化、つまり任意の点で微分がバニッシュしない）な関数φを用いて変換
C(a,b,c)=Σ|φ({y}_{i})-φ(F({t}_{i},a,b,c))|^2

D:所謂1ノルム
D(a,b,c)=Σ|{y}_{i}-F({t}_{i},a,b,c)|

（Q1）このとき、以下の命題のうち、同値な命題はどれとどれですか？

P「(a,b,c)がAの極値点である」
Q「(a,b,c)がBの極値点ある」
R「(a,b,c)がCの極値点ある」
S「(a,b,c)がDの極値点である」

(Q2)重みづけの意味について：
*実際、最急降下法のプログラムを作ってみると、

Dの場合で、直線に近い形状になるように変換した場合（たとえばシクモイドの場合logをφに取る）。
Cの場合で、変化が緩慢なiに重みをつけた場合。
Aの場合。

の順に速度が速く、いずれの場合もだいたいの場合には、まあまあ（10000回ぐらい再起計算すれば）
まあ、見た目に近いグラフが出てきます。

だとしたとき、wやφというのは、何を意味しているのでしょうか？

ここで、最急降下法は、以下の意味で考えている
最急降下法の初期パラメータを(a_0,b_0,c_0）とし、
k回目の計算値を({a}_{k},{b}_{k},c_{k})と記載する。

このとき、({a}_{k},{b}_{k},c_{k})は、以下の漸化式を再帰的に数値計算することで求める。

(a_{k+1},b_{k+1},{c}_{k+1})=(a_{k},b_{k},{c}_{k})-ε*grad(J(a,b,c))

但し。Jは、A,B,Cいずれかの評価関数で、εは充分小さい正定数

(Q4)εのテンソル化：
というほど大げさなものではありませんが、εを正値の対角行列にした場合
収束が早いことがあります。この場合εの異方性がフィッティングパラメータの収束性
どのように作用しているのでしょうか？また、こんなことをしてもいいのでしょうか？

回答No.2

stomachman · ベストアンサー

No.1に付けられたコメントについてです。

> （ア）曲線回帰の問題は、線形化してから回帰すれば、直線回帰と同じなので、線形化してから回帰すればいい。
> なんていうことが、成書にかいてあったりするのですが…。また、線形化してから回帰したほうが
> 収束が早かったりするケースもままあるのですが…。

　収束が早いも何も、線形なら反復なしの一発で答が出るでしょ？

> のε[i]が正規分布となるような回帰のしかたと考えてよろしいでしょうか？

　モデルが線形だろうが非線形だろうが関係なく、
(a) サンプルごとの残差をそれぞれ独立な誤差であるとみなす事ができ、しかも
(b) サンプルごとの誤差の分布が正規分布であって、その分散が既知である
という条件が成立てば、(b)の既知の分散を使って重み付き最小二乗法解を計算すると、それはパラメータの最尤推定になります。

　要するに、ご質問のBが最小二乗法です。ご質問のAはBにおいて重みがすべて1であるという特別な場合に過ぎません。
　ご質問のCは、yではなくてφ(y)こそが測定値であって、それにモデルφ(F)を当てはめるのだと考えれば、Aと同じことです。
　もちろん、Cじゃなくてさらに重みを付けたものを考えるべきです。すなわち、yにFを当てはめる場合の残差 ε = y - F に対して、φ(y)にφ(F)を当てはめる場合の残差 δ = φ(y) - φ(F) は |ε|が小さい時には
　　δ = φ(y) - φ(F) = (φ(F+ε) - φ(F) ) ≒ εφ’(F)　　（ただしφ’(x) = dφ/dx）
となるので、εが平均0分散σ^2の分布に従うならば、δ は平均0分散(σφ’(F))^2の分布に従う。これに応じて重みを付けてやれば良いですね。

　モデルを線形にする話については、ある非線形のモデルFがあるとき、これをうまく変形・変換して線形のモデルGに書き換えられたとします。このとき、Gにおいて(a)(b)が成立てば、Gの最小二乗解がモデルF, G（どっちでも）の最尤推定になっていると言えます。てか、この場合には、Gこそが本来のモデルであって、Fこそがそれをへんてこにいじくったものに過ぎなかったわけです。
　一方、非線形のモデルFにおいて(a)(b)が成立つけれども、Fを変形して作った線形モデルGでは(a)(b)のどちらか、あるいは両方が成立たなくなる、という場合はどうかと言いますと、Gの最小二乗解はFの最小二乗解とは違うから最尤推定ではない。しかしGの最小二乗解はFの最小二乗解にほぼ近いことは確からしい。なので、、Gの最小二乗解を、Fに関する重み付き非線形最小二乗法の繰り返し演算の出発値として利用すれば、（デタラメな出発値に比べて）収束がずっと速い。

> F(t,P,G,E)=P/(1+EXP(G(t-E))

　すなわちパラメータがP,G,Eである非線形モデルFがあって、しかもこのモデルが適切である（(a)(b)が成立つ）という場合を例にして、線形化したモデルを構成してみましょう。
　　y[k] = P/(1+EXP(G(t[k]-E))) + ε[k]
　この場合、サンプル点t[k]が十分沢山あるなら、一度微分方程式に持っていく手があります。y,εをtの連続関数と見れば
　　y’= Gy - GP -εG
であるから、
　　y[k] = G∫y dt - GP∫dt - G∫εdt　　(積分はt[0]～t[k]の定積分）
である。そこで残差を
　　δ[k] = - G∫εdt
として∫y dtを数値積分で計算してやれば線形モデル（パラメータはGと-GP)が得られます。たとえば台形則で数値積分するなら、サンプル点t[k]を小さい順に並べ直しておいて
　　y[k] = (G/2) Σ(t[k]-t[k-1])(y[k]+y[k-1])- GP(t[k]-t[0]) + δ[k]
とやる。すると、δ[k] の２乗和を最小化するG, Pは線形最小二乗法で一発で計算でき、これらが決まると元のモデルからEも決まります。しかし残差は
　　δ[k] = - G∫εdt + (数値積分の誤差)
であるから、「サンプルごとに独立」という条件が満たされませんし、いやそれ以前に、こうして得たP, G, Eはεの二乗和を最小にしないから、Fの最小二乗解ではない。
　けれども、P, G, EはFの非線形最小二乗法の出発値として使える。出発値が解に近ければ、繰り返し計算の途中で数値的不安定性などのおかしなことが起こる恐れも小さい（真の解がそういうヤバい所にある場合は別ですが）。

　つまり、まず適切な出発値を素早くみつけ、それを使って、本来の非線形最小二乗法をきちんと遂行する、ということによって、正しい解を速く見つける、という戦略をとるわけです。

　[話2]以降は非線形最小二乗法の探索の工夫の話でしょうから、(Q4)の話を具体例に適用した個別の工夫ということですね。もちろん、個別のモデルにおいて工夫をするのは結構なことですが、ベースになっているのが最急降下法ではあんまり面白みがありません。
　一方、一般的に利用できる工夫としては：
・Gauss-Newton法において、局所１次近似をして出てくる二階微分をまじめに計算するか、あるいは二階微分は近似式で代用するか、いや微分そのものを差分で代用するか（微分ががとても複雑な式になるモデルでは、繰り返し回数が増えても毎回の計算が速くなるから有用です）。
・繰り返しの度に評価関数が調子良く減少する時にはパラメータの変化を大きく、さもなければ小さくするように自動調節する。
・迷子になったときのために、最もオバカな最急降下法を加味しておく。
・パラメータの空間中で一度調べた所の情報を保存して、周囲の様子を推定する材料にする。
などなど。これらを知った上で個別の工夫を考えると良いでしょう。

お礼

回答ありがとうございました。消化に時間を要しすみません。降下法の部分は、理解が追いついていないので…。まず、「正しいノルム」について整理すると、こんな感じでしょうか？

こんな感じでしょうかと言いながら、自分でも怪しいと思ってますが。

【ノーテーション】

(t[i],y[i])　i=1,2,...,n

を、F(t;a,b,c)でフィッティングする問題を考える。

第i個目のデータの残差を以下のように定義する。

ε[i]=|y[i]-F(t[i],a,b,c)|

となる。

【残差iの確率分布】
正値の関数
w(t):R→R+

によって、

η[i]=\sqrt{w(t[i])}ε[i]

とした場合に、η[i]が、確率分布Pに従うとする。
このとき、

i番目のデータのみに着目した場合、
フィッティングパラメータの値が(a,b,c)である確率p[i](a,b,c)は

p[i](a,b.c)=P(ε[i])=P(|y[i]-F(t[i],a,b,c)|)

となる。

【残差の確率分布】

それぞれの測定が独立試行である場合で、
n個のデータを知っていた場合、
フィッティングパラメータの値が(a,b,c)である確率p(a,b,c)は、

p(a,b,c)={Π}_{i=1}^{n}p[i](a,b,c)

となる。ただし、Πは連乗積である。

【最も尤もらしいa,b,c】
それぞれの測定が独立試行である場合で、n個のデータを知っていた場合、
最も尤もらしいa,b,cとは、p(a,b,c)が最大となるa,b,cである。

【Pが正規分布の場合】

p[i](a,b,c)=N[0,1](η[i])

であるときには、

p(a,b,c)=\frac{1}{{2\pi}^(n/2)exp(-{\sum}_{i=1}^{n} w(t[i])\frac{ε[i]}{2} )}

が、最大となるa,b.cが、最も尤もらしい。

一般に、
N[m,\sigma](x)

は、xが0に近ければ近いほど大きいため、

{\sum}_{i=1}^{n} w(t[i])\frac{ε[i]}{2}

が、最小となる、a,b,cが尤もらしい。

但し、

N[m,\sigma](x)=\frac{1}{\sigma\sqrt{2\pi}}exp(-\frac{(x-m)^2}{2{\sigma}^2})

とする。

回答No.3

stomachman

No.2に付けられたコメントについてです。

　実務上の多くの問題では、誤差の分布がきちんと分かってるのかどうかだなんて、あらまあ、そんな怖い顔してあんまり問いつめないでほしいような状況が多くて、そういう場合に最尤推定にこだわるのは無駄、いやもう間違いだと言っていいと思います。
　けれども、観測精度の限界に近いようなデータに基づいて物を言おうとする場合は、誤差の確率分布の丁寧な考察（t[i]の含む誤差まで含めて）をきちんとやる必要があり、そしてパラメータの値の範囲（危険率の設定によって決まります）をきちんと計算しなくてはならない。

> i番目のデータのみに着目した場合、
> フィッティングパラメータの値が(a,b,c)である確率p[i](a,b,c)

　まず、「フィッティングパラメータの値がある確率分布に従う」という見方は、フィッティングパラメータの計算方法を決めたときにだけ意味のある概念ですね。そこで、具体的に計算方法を（たとえば最小二乗法に）固定して考える必要がある。
　次に、それがどんな計算方法であれ、「全てのサンプルそれぞれにおいて、測定値がランダムな誤差のためにどんな値になるか、ということの影響を受けてフィッティングパラメータが決まる」という、そういう計算方法であるに違いない。だから、簡単に「i番目のデータのみに着目」という訳には行きません。
　計算方法が決めてあれば、「i番目のサンプルのデータy[i]がその計算方法によって）フィッティングパラメータa[j] (j=1,2,...,m)の値にどう影響するか」を関数a[j](y[1],...,y[n])として扱う、ということができ、これがy[i]が従う確率分布とパラメータベクトルa = (a[1],...,a[m])が従う確率分布の関係式に他なりません。
　線形最小二乗法ならa[j]はyの線形結合で書けて、これは正規方程式そのものです。また、それ以外（非線形最小二乗法、あるいは、もっと他の方法）の場合でも、真値z[i]を使ってy[i]=z[i]+e[i] (eは誤差)と書いたとき、パラメータベクトルaをeに関してテイラー展開(a(z+e) = a(z) + (∂a/∂e)e + ...)すれば、|e|が小さいなら線形結合あるいは低次の多項式で近似できることが多いでしょう。

　　ところで、最小二乗法が最尤推定になるかどうかは、独立なサンプルy[i](i=1,2,...)がある場合、yに関する尺度を非線形変換して各サンプルy[i]の誤差の分布が正規分布になるようにできる、そういう変換φが可能かどうか、というところに掛かっているでしょう。それができるなら、
　　Y[i] = φ(y[i]), G(t[i],a,b,c) = φ(F(t[i],a,b,c) )
として、YにFをfittingすれば文句なし。一方、そうは行かない場合には、尤度を残差二乗和で代用するわけには行きませんから、「尤度を最大化するフィッティングパラメータを求む」という極値問題を直に解くしかない。
　こちらもご参考に → http://okwave.jp/qa/q7188529.html

お礼

最適化の部分以外、おおむね理解できたと思います。
懇切丁寧にありがとうございました。

最適化の部分は、再度別の形で質問します。
その折は、よろしくお願いします。

返答が遅れてすみません。

回答No.1

stomachman

残差ε
ε[i] = y[i]-F(t[i],a,b,c)
を「小さく」する、という話でしょう。
　　評価関数A
　　A(a,b,c) = Σ(|ε[i]|^2)
を最小にする(a,b,c)を求めるのは最小二乗法ですね。特に、もしFがa,b,cについて線形なら線形最小二乗法です。
　　B(a,b,c) = Σ(w[i](|ε[i]|^2))　　(w[i]≧0)
を最小にするなら重み付き最小二乗法。
　　δ[i] = φ(y[i])-φ(F(t[i],a,b,c)
　　C(a,b,c)=Σ|δ[i]|^2
を最小にするのも最小二乗法で、データの尺度をφで変換したもの。特に、もしφ(F)がa,b,cについて線形なら線形最小二乗法です。
　これらに比べて
　　D(a,b,c)=Σ|ε[i]|
を最小化するのはかなり毛色が違います。

(Q1)
　残差ε[i]が全て0になる解(a,b,c)が存在するなら命題P～Sは同値ですが、さもなければ命題P～Sは一般に（i.e., 特殊なタマタマを除いて）同値ではないことは明らかでしょう。

(Q2)
　w[i]の意味する所は「サンプルiにおいて残差を相対的にどの程度小さくしたいか」です。w[i]が大きいほど、他のサンプルにおける残差が大きくなるという犠牲を払ってでも、サンプルiでの残差を小さくしたい、ということ。
　サンプルごとにデータyに含まれる誤差が従う確率分布が異なる場合に用います。残差があまり大きくないうちは、w[i]=(分散の平方根の逆数)にするのが最適です。こちらもご参考に→ http://okwave.jp/qa/q3557287.html

　非線形変換φの意味する所は（仰る所の「目的変数」の）尺度の変更です。一つには非線形の問題を線形の問題に（あるいは非線形性の少ない問題に）変換するための手段として、wと組み合わせて使われます。たとえば
　　F(t,a,b,c) = a (t^b) exp(-ct)
で、データyが常に正である場合、φ=ln を使って
　　ln(F(t,a,b,c)) =ln(a) + b ln(t) - ct
のように変換して残差を
　　η[i] = ln(y[i])-ln(F(t,a,b,c))
とします。もしε[i]が従う分布がiに依らないなら、この変換をすると当然、η[i] が従う確率分布がサンプルごとに異なることになるので、これを補正してやるために重みw[i]が必要になります。すなわち
　　η[i]=-ln(F(t,a,b,c)/y[i])=-ln(1-ε[i]/y[i])≒ε[i]/y[i]　(|ε[i]/y[i]| <<1のとき）
であるから、重みw[i] = y[i]^2を使えば、Aの意味での最小二乗法に近い結果を線形最小二乗法で計算できることになります。その他の処方についてはたとえば→ http://okwave.jp/qa/q198558.html
　逆にまた、サンプルごとに誤差が従う確率分布が異なっていて、しかもその分散がサンプルyの値に依って決まる場合に、どのサンプルでも誤差が従う確率分布が同じになるようにするための手段として、非線形変換φを使うこともあります。（この手段を使うのは、最小二乗法に限ったことではありませんが。）→ http://okwave.jp/qa/q5190616.html

> Dの場合で、直線に近い形状になるように変換した場合（たとえばシクモイドの場合logをφに取る）。
　はてさて、Dにはφが入っていませんよ？もしかして、Cの話でしょうか。それにしても、「直線に近い形状」「シクモイドの場合」とは、何と何の関係について仰っている? もしかして、tとyとの関係でしょうか？だとすれば、そんなのφの選択とは無関係です。φは「説明変数」t についてではなく、「目的変数」の変換ですからね。（もちろん、tについて何かの変換をしてもよくて、それはモデルFを変更した、ということに他なりません。）

> Cの場合で、変化が緩慢なiに重みをつけた場合
　　これはBの場合について仰っているのでしょうかね。しかし、「変化が緩慢なi」とは一体どういう意味なのでしょうか。iの番号の付け方はどうであっても構わないわけで、結果にも途中経過にも全く影響がないことは、式を見れば明らかでしょう。

(Q4)εのテンソル化
　（非線形最小二乗法に限らず）非線形の極値問題の数値解法では、要するに極値を探して動き回る。その時に「素直に最急降下方向に動かずに、どっち向きに動くかにバイアスを与える」というのが、おそらく「εのテンソル化」と仰る意味でしょう。
　最急降下法は「遅くてもいいから簡単な計算式・プログラムで済ませたい」という時に使われます。しかし、「（最急降下方向の情報も場合によっては利用するけど）もっと他の情報を利用して、動く方向と量を決める」というもうちょっと知恵を使った手法がいろいろあります。たとえばMarquardt法等を参照なさると参考になるかも。

　ところで、（ご存知かも知れませんが）非線形極値問題では、出発値の選び方が悪いとlocal minimimに落ち込んで誤った結果に収束してしまうことはしばしばあります。（最急降下法を含め）ちょこっとずつ動いてみるというやり方では（局所しか見ていないので）この問題は避けられません。これを（なるべく）避けるには、上記φによる変換などの手法で問題を線形化して適切な出発値を計算するか、あるいは「焼き鈍し法」などのlocal minimumを抜け出すメカニズムを入れる必要があります。

お礼

返信が遅くなりすみません。
丁寧なご回答ありがとうございます。

回答の理解がおいつかず、今もまだ全部は理解しきれていませんが、理解した？とこまで書きます。

(Q1)
＞残差ε[i]が全て0になる解(a,b,c)が存在するなら命題P～Sは同値ですが、さもなければ命題P～Sは一般に（i.e., 特殊なタマタマ＞を除いて）同値ではないことは明らかでしょう。

ですよね。という言い方は失礼かもしれませんが。

かなり毛色が違うＤ（命題Ｓ）はともかくとして、（質問にあげておいてなんですが、一旦Ｄについては、忘れます。）

よく、

（ア）曲線回帰の問題は、線形化してから回帰すれば、直線回帰と同じなので、線形化してから回帰すればいい。

なんていうことが、成書にかいてあったりするのですが…。また、線形化してから回帰したほうが
収束が早かったりするケースもままあるのですが…。

ノイズレベルなりなんなりについて考察したうえで線形化したほうがよいという結論が
先にあってしかるべきで、その結果として、「正しい」回帰の仕方が決まる
（A～Cのいずれかに定まる）という理解でよいでしょうか？

また、若干話が飛びますが、「正しい回帰の仕方というのは、

【Q5】（新たな質問）

のε[i]が正規分布となるような回帰のしかたと考えてよろしいでしょうか？

(Q2)
>サンプルごとにデータyに含まれる誤差が従う確率分布が異なる場合に用います。残差があまり大きくないうちは、w[i]=(分散の平方>根の逆数)にするのが最適です。こちらもご参考に→ http://oshiete.goo.ne.jp/qa/3557287.html

ε[i] = y[i]-F(t[i],a,b,c)
B(a,b,c) = Σ(w[i](|ε[i]|^2))　　(w[i]≧0)

のε[i] が、いずれのiにおいても、正規分布に従い、かつ、σ[i]がiごとに違う場合には、

χ^2 = Σ([y[i]-F(t[i]; a, b, ...)]/σ[i])^2

を評価関数として取るのが公平であるということでしょうか？
【ご参考】

http://oshiete.sponichi.co.jp/qa5841569.html

のANo.2のhitokotonusi 氏の回答。

>もしε[i]が従う分布がiに依らないなら、この変換をすると当然、η[i] が従う確率分布がサンプルごとに異なることになるので、これを補正してやるために重みw[i]が必要になります。

η[i]=φ({y}_{i})-φ(F({t}_{i},a,b,c))

誤差が、ε[i]が従う分布が正規分布ですらない場合には、
φは、η[i]が正規分布となるようにφを選ぶのが望ましい
ということでしょうか？

また、

η[i]=が、従う確率分布がいずれも正規分布である場合で、標準偏差σ[i]も等しい場合
⇒新たに重みをつける必要がない。

η[i]=-ln(F(t,a,b,c)/y[i])が、従う確率分布がいずれも正規分布である場合で、標準偏差σ[i]も等しくない場合
⇒Σ(η[i]/σ[i])^2

を最小化するのが公平であるということでしょうか？

＞はてさて、Dにはφが入っていませんよ？もしかして、Cの話でしょうか。

すみません。私のミスです。

それにしても、「直線に近い形状」「シクモイドの場合」とは、何と何の関係について仰っている?

シクモイドは、

F(t,P,G,E)=P/(1+EXP(G(t-E)))

を考えています。

>もしかして、tとyとの関係でしょうか？

すみません。いろいろな話をごちゃごちゃにしてしまいました。以下の３つの話にわけます。

[話1]φ=LNとして　モデルCに基づいた回帰をすると、収束が早い。
[話2]以下の漸化式に
(a_{k+1},b_{k+1},{c}_{k+1})=(a_{k},b_{k},{c}_{k})-c*grad(J(P,G,E))

に対し、Ψ(P,G,E):R^3→R^3 を用いて、
(a_{k+1},b_{k+1},{c}_{k+1})=(a_{k},b_{k},{c}_{k})-c*grad(Ψ(P,G,E))

とする。但し。Jは、A,B,Cいずれかの評価関数で、cは充分小さい正定数

[話3]変化が緩慢なiに重みをつけた場合 (後述)

話1は、だいたいわかったような気がします。
話2は、特に、

(a_{k+1},b_{k+1},{c}_{k+1})=(a_{k},b_{k},{c}_{k})-c*grad(|LN(P)|,G,E)

とすると、収束が早くなるように思います。

> Cの場合で、変化が緩慢なiに重みをつけた場合
　　これはBの場合について仰っているのでしょうかね。しかし、「変化が緩慢なi」とは一体どういう意味なのでしょうか。

シクモイドの場合は、「シクモイドの二階微分の極大値を与えるt」よりも
小さいtについて、重みをつけると収束が早いように思います。

例えば、

(100/(1+EXP(-(t-10))))

が「最適」のケースでは、
t<9

のデータに重みをつけると、収束が早いように思います。

(Q4)εのテンソル化
＞（非線形最小二乗法に限らず）非線形の極値問題の数値解法では、要するに極値を探して動き回る。その時に「素直に最急降下方向＞に動かずに、どっち向きに動くかにバイアスを与える」というのが、おそらく「εのテンソル化」と仰る意味でしょう。

はいそうです。

＞　最急降下法は「遅くてもいいから簡単な計算式・プログラムで済ませたい」という時に使われます。しかし、「（最急降下方向の＞情報も場合によっては利用するけど）もっと他の情報を利用して、動く方向と量を決める」というもうちょっと知恵を使った手法が＞いろいろあります。たとえばMarquardt法等を参照なさると参考になるかも。

一次近似ではなく２次近似をつかうといった感じでしょうか？

二次近似はいいのですが、
一次近似なのに、異方性を持たせたら収束がはやくなるとか、
上の[話2]のケースとかがなんか気持ちが悪いのです。

(a_{k+1},b_{k+1},{c}_{k+1})=(a_{k},b_{k},{c}_{k})-α*grad(|LN(P)|,G,E)

の極限として定まる微分方程式が、同じ極値（最適とは限らなくてよい）に漸近するかというのが、
気になるといいますか…。