データ列(xi,yi)(i=1,...,n)を関数
y=α(x+γ)^2 exp(-β(x+γ)) (α>0,β>0,γ>0)
でフィッティングしたいです。対数をとって普通に最小2乗法で解こうとしたら得られる連立方程式が線形でなくて解けませんでした。
どうしたらいいのでしょうか?

このQ&Aに関連する最新のQ&A

A 回答 (3件)

完全な意味での最小二乗法ではないけれど、無理矢理線形の問題に帰着して簡単に計算するstomachman我流のやりかたがあります。

(超平面法と勝手に呼んでいます。)この方法は、「データがモデルに良く合う」という時に限り「最小二乗解とほぼ一致する」という性質があります。

モデルをデータと区別して、
f(t) = α((t+γ)^2) exp(-β(t+γ))
と書きましょう。
両辺をtで微分すると
f'(t) = 2α(t+γ) exp(-β(t+γ))-βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2α(t+γ)^2) exp(-β(t+γ))-(t+γ)βα((t+γ)^2) exp(-β(t+γ))
(t+γ)f'(t) = 2f(t)-(t+γ)βf(t)
(t+γ)f'(t) = (2-βγ)f(t)-βtf(t)
という微分方程式が得られます。
 この両辺をtで積分すると、(t=x1~xiまで)
∫(t+γ)f'(t)dt = (2-βγ)∫f(t)dt-β∫tf(t)dt
左辺を部分積分して
(xi+γ)f(xi)-(x1+γ)f(x1)-∫f(t)dt = (2-βγ)∫f(t)dt-β∫tf(t)dt
整理すると
xif(xi)-x1f(x1) = (3-βγ)∫f(t)dt-β∫tf(t)dt-γ(f(xi)-f(x1))
です。
ここで、f(xi)の代わりにyiを使っても、「データとモデルが良く合う」のなら大差は出ない筈。だから
∫f(t)dtの代わりにyをx1~xiまで数値積分したものJi、
∫tf(t)dtの代わりにxyをx1~xiまで数値積分したものKi、を用いることにして、残差をεiと書くことにし、
xiyi-x1y1 = (3-βγ)Ji-βKi-γ(yi-y1)+εi (i=1,2,....,n)

改めて
Yi = xiyi-x1y1
Li = yi - y1
A = (3-βγ)
B = -β
C = -γ
とおくと、
Yi = A Ji + B Ki + C Li + εi (i=1,2,....,n)
という線形問題になったわけで、これは簡単に解けます。
 ところで、A,B,Cはβとγだけで表される一方、αが出てきませんね。
つまり
・A,B,Cからβ,γを求めると1通りに決まりません。これは気にしないで、たとえばAを無視してβ,γを計算します。ついでに、こうして得たβ,γを使って(3-βγ)を計算してみます。もしこれがAとまるで合わないようなら、「データがモデルに良く合う」という条件が満たされていない。yiに含まれるノイズが大きすぎるんです。
・αを求めるには、
f(t) = α((t+γ)^2) exp(-β(t+γ))
を再び使います。γとβは既に分かっているから、簡単な線形最小二乗法の問題ですね。

より高精度の最小二乗解を求めるには、超平面法で得た解を初期値にして、非線形最小二乗法(ガウス・ニュートン法など)を使って改良すると良いです。(教科書としてはいつも、中川・小柳「最小二乗法による実験データ解析」東京大学出版会をお薦めしています。)
    • good
    • 0

最小二乗法と非線型2情報の話しは既に出ていますので.ORの話しにしましょう。



最適化問題で非線型問題がどうしても融けない場合に使います。
まず.でたらめにアルファ・べーた・ガンマーの値を決めます。
簡単にするのでしたらば対数で1e-8から1え+8の間が均等になるような点(1.10.100とか標準数(JISZ番号忘れ)とかあります)を選び.残差を求めます。
そのなかで.もっとも残差が少ない.数値の組み合わせを選びます。

次に.どれか一つをちょっと変えて計算してみます。「ちょっと」の程度ですが.切りの良い(2進数で考えてください)1/128ぐらいではじめます。もし.最初の計算よりも残差が少なくなったらばその解をえらびます。もし.残差がおおきくなってしまったらは.「ちょっと」の程度を更に1/2にします。計算していて.計算機Eの問題が出てきたら計算を止めます。
この考え方を基にして.3つの変数の解の付近(目安としては前回の計算した範囲程度)を128分割して演算をしてその中のもっとも残差の少ない点を解とします。

最低でも倍精度演算が必須です。誤差関数を求めるために8倍精度とか無限倍精度(文字演算で行うために計算がやたら遅いのでやりたくない)とかを使う方法もあります。

OR問題の亜種です。私は場当たり法と呼んでいます。メッシュで計算する法ですから.解の巣問題や鞍部問題(1の方の参考文献を読んでください)に対して(計算を間違えなければ)強力な計算方法です。しかし.計算時間がやたらかかるので.評判はよろしくないです。8-10変数で.100点.80486で8時間なんてざらです(真の最適解とは限らないので.初期値をちょっとずらして再計算するとすぐに1-2ヶ月使ってしまいます)から。
    • good
    • 0

これは非線形最小二乗法の問題になりますが,


過去いくつかのスレッドがあり,処方もいろいろ回答されています.
たとえば,
http://oshiete1.goo.ne.jp/kotaeru.php3?q=179919
http://oshiete1.goo.ne.jp/kotaeru.php3?q=190632
http://oshiete1.goo.ne.jp/kotaeru.php3?q=98446
http://oshiete1.goo.ne.jp/kotaeru.php3?q=97271
http://oshiete1.goo.ne.jp/kotaeru.php3?q=33318
など.
他にも質問検索で「最小二乗法」「最小自乗法」とやるとかなりヒットします.
    • good
    • 0

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Q最小二乗法の過程で分散共分散行列が・・・・

最小二乗法の過程で分散共分散行列が目にしました。

重回帰モデルが
y=Xβ+ε
で与えられるとき、最小二乗法を施すと
βの期待値bが
(XT・X)-1・XT・yで表せるのですが、
(ただし、XTは行列Xの転置行列、-1は逆行列を表します。)

大事なのがここからで、
この (XT・X)-1  というものが、有名な形らしく、分散共分散行列と呼ばれるらしいのです。
どうして、この行列が分散を表す行列になるのかが、いまいちつかめないのです。ご存知の方がいらっしゃいましたらぜひ教えてください。
助けてください!!よろしくお願いします!!

Aベストアンサー

#1です。再度補足。

1)分散の定義そのものです。
詳しくは確率変数の平均値が0でないときの分散の定義を見てください。

2)(X'X)^(-1) X'Xβ- β=0になるということでしょうか?
その通りです。(X'X)^(-1) X'X = I 、すなわち単位行列でしょう?

3)E[(b-β)(b-β)'] = σ^2 (X'X)^(-1)
地道に計算すれば解けるのですが、
E[(b-β)(b-β)'] = E[((X'X)^(-1) X'(Xβ + ε) - β)((X'X)^(-1) X'(Xβ + ε) - β)']
= E((X'X)^(-1) X'ε) ((X'X)^(-1) X'ε)']
= E[((X'X)^(-1) X'ε) (ε'X((X'X)^(-1))]
= ((X'X)^(-1) X'E[εε']X((X'X)^(-1))
= ((X'X)^(-1) X'σ^2 X((X'X)^(-1)) ←σ^2 はスカラーなので前に出せる(単位行列は消す)
= σ^2 (X'X)^(-1) X'X(X'X)^(-1)
= σ^2 (X'X)^(-1)
となります。

#1です。再度補足。

1)分散の定義そのものです。
詳しくは確率変数の平均値が0でないときの分散の定義を見てください。

2)(X'X)^(-1) X'Xβ- β=0になるということでしょうか?
その通りです。(X'X)^(-1) X'X = I 、すなわち単位行列でしょう?

3)E[(b-β)(b-β)'] = σ^2 (X'X)^(-1)
地道に計算すれば解けるのですが、
E[(b-β)(b-β)'] = E[((X'X)^(-1) X'(Xβ + ε) - β)((X'X)^(-1) X'(Xβ + ε) - β)']
= E((X'X)^(-1) X'ε) ((X'X)^(-1) X'ε)']
= E[((X'X)^(-1) X'ε) (ε'X((X'X)^(-1))]
= ...続きを読む

Q異なる4点A(α)、B(β)、C(γ)、D(δ)で、|α|=|β|=|γ|=|δ|、α+β+γ+δ=

異なる4点A(α)、B(β)、C(γ)、D(δ)で、|α|=|β|=|γ|=|δ|、α+β+γ+δ=0のとき、A、B、C、Dを頂点とする四角形が長方形になることの証明を、どなたかお願いします。

Aベストアンサー

(1) 2次元ユークリッド平面上のベクトルの話だという限定を付けないと、長方形にはならない。(3次元なら、たとえば原点に重心がある正四面体の頂点がα,β,γ,δでも条件を満たすでしょ。)
(2) |α|=0の場合は例外だし、α,β,γ,δのうちに同じものが含まれる場合も例外。
ということに注意した上で
(3) |α|=|β|=|γ|=|δ|=1の場合に証明すれば、他の場合は自明なので、=1の場合だけ考える。
(4) x = (α+β) とすると、αとxがなす角θはxとβがなす角と同じ。
(5) (γ+δ) = -xでなくちゃならない。で、γとxがなす角ξはxとδがなす角と同じ。
あとはθ=ξを示せばよかろ。

Q最小二乗法 擬似逆行列

下のサイトの説明を読んで最小二乗法の勉強をしています。

未知パラメータxが1個(N=1)、出力yが10個(M=10)のとき、
擬似逆行列を求めようとすると、(A#A)^-1 が1×1行列になってしまいます。
このとき逆行列はどのように求めたら良いのでしょうか??勉強不足ですみませんが、よろしくお願いします。

http://www.star.t.u-tokyo.ac.jp/~kaji/leastsquare/leastsquare_main.htm

Aベストアンサー

転置を'で表すことにします。

1×1行列 A'A = [a] の逆行列は [1/a] ですので、擬似逆行列は

(A'A)^(-1)A'=[1/a] A'=(1/a) A'

と計算されます。蛇足ながらこの場合の擬似逆行列は1×M行列となります。

Qax^3+bx^2+cx+d=0がα、β、γを解に持つならばax^3+

ax^3+bx^2+cx+d=0がα、β、γを解に持つならばax^3+bx^2+cx+d=a(x-α)(x-β)(x-γ)と変形できることを示せ。

Aベストアンサー

因数定理を証明しろってことなのかな?

z の多項式 f(z) が根 z=0 を持つとすると、
f(z) は 定数項が 0。よって、z で割り切れる。
定数項が 0 であることは、
f(z) を z の降冪または昇冪に整理して、
z=0 を代入してみれば判る。

x の多項式 F(x) が根 x=α を持つ場合は、
F(α+z) を z の多項式と見て、
上の補題を適用すれば解る。
F(α+z) が根 z=0 を持ち、z で割り切れるので、
z = x-α を代入すれば、
F(x) は x-α で割り切れている。

F(x) が x-α で割り切れれば、
F(x) / (x-α) は x の多項式である。
この多項式が根 x=β を持つとき…

…と繰り返してゆけば、
n 次多項式は n 個の一次式で割り切れ、
最後の商は定数式になる。
(根の存在自体は代数学の基本定理によるが、
質問の例では解の存在が仮定されているから、
その点は気にしなくてよい。)

最後の商を何か未定係数で置いて
一次式の積を展開してみれば、
最高次の係数の比較から、それが a であると判る。

証明の流れを見れば解るように、
これは、α,β,γ の中に同じものがあっても、
それを重解とみなせば、成り立つ。

因数定理を証明しろってことなのかな?

z の多項式 f(z) が根 z=0 を持つとすると、
f(z) は 定数項が 0。よって、z で割り切れる。
定数項が 0 であることは、
f(z) を z の降冪または昇冪に整理して、
z=0 を代入してみれば判る。

x の多項式 F(x) が根 x=α を持つ場合は、
F(α+z) を z の多項式と見て、
上の補題を適用すれば解る。
F(α+z) が根 z=0 を持ち、z で割り切れるので、
z = x-α を代入すれば、
F(x) は x-α で割り切れている。

F(x) が x-α で割り切れれば、
F(x) / (x-α) は x の多項式である。...続きを読む

Q最小二乗法の過程で分散共分散行列が・・・

最小二乗法の過程で分散共分散行列が目にしました。

重回帰モデルが
y=Xβ+ε
で与えられるとき、最小二乗法を施すと
βの期待値bが
(X'・X)-1・X'・yで表せるのですが、
(ただし、X'は行列Xの転置行列、-1は逆行列を表します。)

このbの分散共分散行列を計算すると
E[(b-E(b))(b-E(b))'] = σ^2 (X'X)^(-1)
となるのですが、
この計算をする際にE(b)=β
とするのですが、どうしてこれが成り立つのでしょうか?
教えてください!!よろしくお願いします。

また、もしよろしければ、次のことを教えてください。ご存じなければ、上の質問だけでももちろん助かります。

bの分散共分散行列
E[(b-E(b))(b-E(b))']
を計算すると
E[εε']=σ^2I
となる過程が出てくるはずです。(Iは単位行列)
どうしてこの変形が出来るのでしょうか?
単位行列にはならないはずですが、単位行列になるという
仮定を用いているようです。
どうしてこの仮定が生じたのか、ご存知の方いらっしゃれば教えてください。とても困っています。

最小二乗法の過程で分散共分散行列が目にしました。

重回帰モデルが
y=Xβ+ε
で与えられるとき、最小二乗法を施すと
βの期待値bが
(X'・X)-1・X'・yで表せるのですが、
(ただし、X'は行列Xの転置行列、-1は逆行列を表します。)

このbの分散共分散行列を計算すると
E[(b-E(b))(b-E(b))'] = σ^2 (X'X)^(-1)
となるのですが、
この計算をする際にE(b)=β
とするのですが、どうしてこれが成り立つのでしょうか?
教えてください!!よろしくお願いします。

また、もしよろしければ、次...続きを読む

Aベストアンサー

> (X'・X)-1・X'・y
E[(X'X)^(-1) X'y] = E[(X'X)^(-1) X'(Xβ + ε)]
= (X'X)^(-1) X'Xβ + (X'X)^(-1) X'E[ε]
= β + 0 = β


最小自乗法でおかれる通常の仮定の中に
E[ε] = 0
E[εε'] = σ^2 I
が含まれています。この仮定をゆるめたモデルも存在しますが、最初はここが基本です。したがってこの仮定から出てきます。

ちょっとだけ解説すると、前者はモデルが正しい (y = Xβ + ε が期待値の意味で成り立つ) ことを主張しており、後者は誤差項はそれぞれ独立で同じ分散を持つ、ということを主張しています。

Qa+b+c=2で、a>0,b>0,C>0

a+b+c=2で、a>0,b>0,C>0

のときに、a^3+b^3+c^3の最小値を出せ

という問題ってどうやってときますか?

僕が考えたのが、c=2-(a+b)を代入して、aとbそれぞれで平方完成を考えたのですが、式が複雑になります。スマートに解く方法てあるのですか?

Aベストアンサー

コーシーシュワルツの不等式を使うのは、禁じ手でしょうか。

Q最小二乗法 行列

現在以下のページを参考に最小二乗法の勉強をしています。
誤差の二乗ノルムを求めるときに、
 || e || ^2 = e*e
= (y-Ax)*(y-Ax)
= (y*-x*A*)(y-Ax)
= y*y - y*Ax -x*A*y + x*A*Ax  (1)
この次に
       = y*y - 2y*Ax + x*A*Ax
と変形できるのはなぜなんでしょうか?
(1)の第3項目 x*A*y が第2項目と等しくなる過程が分かりません。
後、次の二乗ノルムを微分する過程も良くわかりません。

すみませんが、よろしくお願いします。


http://www.star.t.u-tokyo.ac.jp/~kaji/leastsquare/leastsquare_main.htm

Aベストアンサー

共役転置行列の定義を確認し、公式
(AB)* = (B*)(A*),
A** = A
を証明してみてください。
これにより、(x*)Ay の共役転置は、
{ (x*)(A*)y }* = (y*)(A**)(x**) = (y*)Ax です。

(x*)Ay は、スカラー(1行1列の行列)ですから
{ (x*)Ay }' = (x*)Ay となりますが、
実数値であれば { (x*)Ay }~ = (x*)Ay なので、
結局、
{ (x*)Ay }* = { (x*)Ay }'~ = { (x*)Ay }~ = (x*)Ay です。

Qx+y+z=0,2x^2+2y^2-z^2=0のとき,x=yであることを証明せよ。

クリックありがとうございます(∩´∀`)∩

 ★x+y+z=0,2x^2+2y^2-z^2=0のとき,x=yであることを証明せよ。

この問題について説明をお願いします。

Aベストアンサー

おおざっぱな説明になりますが、左の式を
z=-x-y
として、それを右の式のzに代入します。
それを展開してまとめると
x^2-2xy+y^2=0
という式になります。
あとはこれを因数分解すれば
(x-y)^2=0
となるので、x=yという答えがでます。
与えられた条件がほかになければこれでいいはずです。

Q分子構造の重ねあわせ (行列、最小二乗法など?)

すいません、化学っぽい内容の質問なのですが、
アルゴリズムは数学になると思うのでこちらで質問させていただきます。
また、説明不足を補うために一部プログラムっぽい表記をさせていただきますのでご了承ください。

現在、分子構造を比較するためのプログラムを作成しています。
そのプログラムに構造の「重ね合わせ」機能を追加したいのですが
そのアルゴリズムに最小二乗法などが必要らしく、悩んでいます。
(詳細は長くなるので下に書きます)

分かる方がいましたらご教示お願いします。
よろしくお願いします。


■重ね合わせについて
分子構造のデータは、各原子がそれぞれX, Y, Z座標の3つの値を持っています。
例えば、200原子からなる分子のデータ構造は以下のようになります。

Atom { double x; double y; double z; }
Atom atom = new Atom[200];

重ね合わせる2つの構造は原子数は同じとします。なので、分子構造2つをa1, a2とし、
対応する点の誤差をとっていくと、誤差の平均を求める計算は以下のようになります。

Atom a1 = new Atom[200];
Atom a2 = new Atom[200];
double d = 0;
for (i = 0; i < a1.length; i++)
d += sqrt(pow((a1[i].x-a2[i].x),2) + pow((a1[i].y-a2[i].y),2) + pow((a1[i].z-a2[i].z),2));
d = d / a1.length;

このdの値が最小値をとるようにa2の原子の座標を変更することを「重ね合わせ」と定義します。


■重ね合わせのアルゴリズムについて
基準となる構造に対し、もう一つの構造を並進や回転によって合わせる感じになると思います。
今回質問する箇所はこのアルゴリズムになります。
この際に行列計算や最小二乗法を使うことになるかと思われます。

■質問者の考え
途中まで、しかも間違っている可能性がありますが私の考えたアルゴリズムを追記しておきます。
(i) a1, a2の重心を求める
double x1 = y1 = z1 = x2 = y2 = z2 = 0;
for (i = 0; i < a1.length; i++) {
x1 += a1[i].x; y1 += a1[i].y; z1 += a1[i].z;
x2 += a2[i].x; y2 += a2[i].y; z2 += a2[i].z;
}
x1 = x1 / a1.length; y1 = y1 / a1.length; z1 = z1 / a1.length;
x2 = x2 / a1.length; y2 = y2 / a1.length; z2 = z2 / a1.length;

a1の重心(x1, y1, z1)とa2の重心(x2, y2, z2)が求まります。

(ii) 重心の誤差を求め、誤差をa2の各原子の座標に反映させる
xd = (x2 - x1); yd = (y2 - y1); zd = (z2 - z1);
for (i = 0; i < a1.length; i++) {
a2[i].x -= xd;
a2[i].y -= yd;
a2[i].z -= zd;
}

この方法で誤差はかなり小さくなります。しかし、最小値にはなってないと思われます。
なので、この方法が正しいのかよく分かりません。

すいません、化学っぽい内容の質問なのですが、
アルゴリズムは数学になると思うのでこちらで質問させていただきます。
また、説明不足を補うために一部プログラムっぽい表記をさせていただきますのでご了承ください。

現在、分子構造を比較するためのプログラムを作成しています。
そのプログラムに構造の「重ね合わせ」機能を追加したいのですが
そのアルゴリズムに最小二乗法などが必要らしく、悩んでいます。
(詳細は長くなるので下に書きます)

分かる方がいましたらご教示お願いします。
よろし...続きを読む

Aベストアンサー

こんにちは.
ソースコードは良く見てませんが…
要するに合同であることが期待されるn点のデータ集合XとX'に関して,
Σ ||Xi' - (R*Xi+T) ||^2 -> min
なる回転と並進ベクタを計算したいと考えます.
ここで,XiとXi'はi番目の座標を意味する3次元ベクタであり,
総和演算は範囲(i=1,...,n)とします.

質問者様の方法では上記の変換(R,T)に関して,
並進量Tまでしか解決できません.
ここからさらに回転量Rを解決する必要があります.
基本的な戦略は相関行列を特異値分解してRを決定します.
Rの推定には質問者様の方法で重心を合わせたXとX'に関して
M = Σ Xi'*Xi^T
なる3*3相関行列を計算し,M -> U*D*V^T と特異値分解します(^Tは転置記号).
最小二乗の意味で最適なRは
R = U*V^T
で計算できます.
ただし,回転の向きは適宜調整してください.

データ集合XとX'が相似の意味で同じである場合は
重心を会わせた後にスケール調整を行ってください.

こんにちは.
ソースコードは良く見てませんが…
要するに合同であることが期待されるn点のデータ集合XとX'に関して,
Σ ||Xi' - (R*Xi+T) ||^2 -> min
なる回転と並進ベクタを計算したいと考えます.
ここで,XiとXi'はi番目の座標を意味する3次元ベクタであり,
総和演算は範囲(i=1,...,n)とします.

質問者様の方法では上記の変換(R,T)に関して,
並進量Tまでしか解決できません.
ここからさらに回転量Rを解決する必要があります.
基本的な戦略は相関行列を特異値分解してRを決定します.
Rの推...続きを読む

Q二次方程式x2-3x+4=0の2つの解をα、βとするちき、α+1,β+

二次方程式x2-3x+4=0の2つの解をα、βとするちき、α+1,β+1を解とすると二次方程式を作れ 


複素数の範囲で考えて、次の方程式を因数分解しなさい
x2-x+3

x>yならば、6x+2y>3x+5yであること証明せよ
        

a>0のとき、a+a分の1>2を証明せよ


これらの問題がどうやっても解けません 特に証明系が最も難しいです 何かコツなどがありましたら教えてもらないでしょか? すみません おねがいします

Aベストアンサー

ちなみに1問目で左辺の正式をf(x)とおくと
f(x) = 0の解がαとβだから
f(α) = 0,f(β) = 0が成り立ちます。

なのでα+1, β+1を解にしたかったら
f(x)をf(x - 1)とすればいいです。

なぜならxにα+1, β+1を代入すれば
f(α) = 0,f(β) = 0が成り立つからです。


人気Q&Aランキング

おすすめ情報