Feb 21, 2009

ラグランジュの未定乗数法

統計学を復習していて、どうやらそれを使わないと非常に難しいらしい問題に突き当たったので、ラグランジュの未定乗数法を復習する。
その響きで一瞬怯んでしまうが、前提条件にさえ気を付ければ使う分には何も難しくない、ロピタルの定理のような便利技である。

どこにでも書いてあるので、わざわざここに書くまでもないが、ラグランジュの未定乗数法とは、x1...xnの関数fが束縛条件g=0の下で極値を持つ時、その極点のx1...xnは、f-λg(λは仮置きの定数)をx1...xnで偏微分したものが0になる値として求めることができるというものである。
数式で書くと、

である。ロピタルの定理同様、考えてみるとミラクルである。

例題として、2008年2月1日にフジテレビから放送された「たけしのコマネチ大学数学科」の第75回の問題を取り上げる。

問:

このような円柱と三角錐を組み合わせた鉛筆形の立体があり、体積が決まると、円柱の半径と円柱・円錐の高さは鉛筆形の表面積が最小になるように決まるとする。円錐部分の稜線の長さ(頂点から円錐面に沿って引いた線の長さ)が3の時、円柱部分の半径はいくらか。

解答:
円柱の半径をx、円柱の高さをy、円錐の稜線の長さをzとすると、表面積Sと体積Vは、中学校で習った公式を思い出すと
S=2¥pi x^2+2¥pi xy+¥pi xz
V=¥pi x^2y+¥frac{¥pi x^2}{3}¥sqrt{z^2-x^2}
と求まる。

まずは失敗の例から。普通に考えると、例えばVの式をyについて解いて、Sの式に代入してyを消して、
S=¥frac{2¥left( 3V-¥pi {x}^{2}¥sqrt{{z}^{2}-{x}^{2}}¥right) }{3x}+¥pi xz+¥pi {x}^{2}
このSを最小にするx,zを求めれば良いのだが、よく見ると1つの式に対して変数がx,z,Vと3つある上、例えばxについて偏微分すると
¥frac{¥partial S}{¥partial x}=-¥frac{¥sqrt{{z}^{2}-{x}^{2}}¥,¥left( 6¥,V-3¥,¥pi ¥,{x}^{2}¥,z-6¥,¥pi ¥,{x}^{3}¥right) +2¥,¥pi ¥,{x}^{2}¥,{z}^{2}-4¥,¥pi ¥,{x}^{4}}{3¥,{x}^{2}¥,¥sqrt{{z}^{2}-{x}^{2}}}
(右辺はMaximaの出力)となり、いわゆる強靭な計算力が必要な状況に陥る。ちなみに、Maximaでz,Vについても偏微分させてそれぞれ=0の連立方程式を解かせると、"system too complicated."というエラーを出してgive upする。

そこでラグランジュの未定乗数法である。
Vという条件の下でSの最小値を求める問題なので、

である。(正確にはVは(V-体積)だが、偏微分で体積が消えるので不要)
yについての式を実際に偏微分すると

となり、λ=2/xが得られる。次に、zについての式を偏微分してλを消去すると、

となり、これにz=3を代入すると、x=√5と求まる。
Vもyもわからないままで良い。なんということであろうか。

See more ...

Posted at 22:17 in 数学 | WriteBacks (5)
WriteBacks

Maximaの使い方を探していて、立ち寄りました。 ラクランジュ未定乗数法は確かに簡単ですね。 ただ、今回の場合は、最初の方法では解き方を間違っているために、おかしなことになっています。 単に、xを求めるだけであれば、Vが外生変数ですので、問題を解く時には定数で、 あるVに対するx、zの関係がでてきます。Vで偏微分して、解こうとしたら解けません。 1行目から4行目までは基本的に一緒です。 SYというのは、体積をYについて解いたものを代入した式です。 ここで、体積がzに対して最小化するように決まるのですから、zで偏微分してxについて解けば x=5^(1/2)/3zが求まりますので、z=3をいれれば、x=5^(1/2)となります。 S: %pi*x^2 + 2*%pi*x*y + %pi*x*z; solve(%pi*x^2*y + %pi*x^2*sqrt(z^2-x^2)/3=V,y); SY:%pi*x*z+2*%pi*x*(3*V-%pi*x^2*sqrt(z^2-x^2))/(3*%pi*x^2)+%pi*x^2; assume(z>0); solve(diff(SY,z)=0,x);

Posted by ぼんど at 01/31/2010 05:03:07 PM

ぼんど様ご指摘ありがとうございます。ほんとですね。Vはxやzと独立ではないので、Vで偏微分するのは意味が無いし、ここではVは定数扱いですね。紙の上でyを消したSの式(SYに相当)の形を見て計算不能と決めつけてましたが、これをzで偏微分するとVが消えるので、確かにdiff(SY,z)=0として解けますね。しかし、たまたまそうであっただけで、もしSYがxやzで偏微分してVが消えない形なら、やはりSYを最小にするx,zを求めるのは困難なのでは…

Posted by ynomura at 01/31/2010 10:18:39 PM

Vは与えられる定数で、その中で解きますので、Vが残った時は、V=の式を作って、それを連立させて消すことで、xとzとの式に帰着できます。 dS/dx=(d/dx) f(x,z,V)=0、dS/dz=(d/dx) f(x,z,V)=0となり、どちらかの式からV=を作り、これを使ってVを消去して、x=g(z)という式を導けばOKです。 ちなみに、未定乗数法というのは、代入法を簡単にしたと考えることができます。 未定乗数法が真価を発揮するのは、与えられた条件が陽に解けない場合です。 f(x,y,z)=0という式があったときに、これを、y=f(x、z)といった形に変形することを陽といいます。f(x,y,z)=0は、各変数の関係が式の中に隠れていますので陰(関数)です。 今回の問題の場合は、y=の形に式が変形できていますので、未定乗数法を使わなくても解くことが容易です。 ところが、式によっては、y=という形に式を変形して、代入できない場合があります。 こういうときは、未定乗数法を使うことで、y=という式を意識しないで問題を解くことができます。(ただし、陰関数定理によって、具体的には陽にできないけれども、理論的にはy=が存在することは保証されないといけません) 例えば、z=f(x,y)という式の極値を、g(x,y)=0という条件下で考える時、条件を書き換えて、y=j(x)という形にできるなら、これを代入することでz=f(x,j(x))という関数になりますので、これを微分して dz/dx=df/dx+(df/dy)(dj/dx)=0・・・(1)   とすれば、xが求まります。 ここで、面白いのは、y=j(x)としなくても、g(x,y)=0を用いて未定乗数法を行うと、dj/dxがj(x)という関数を使わずに表現することができます。 j(x)がどんな形かわわかりませんが、あるという前提で、g(x,j(x))=0とします。これを微分すると、dg/dx+(dg/dy)(dj/dx)=0となりますので、式変形すると dj/dx=-(dg/dx)/(dg/dy)です。 この結果を、(1)に代入すると dz/dx=df/dx+(df/dy)(dj/dx)=df/dx-(df/dy){(dg/dx)/(dg/dy)}=0 となります。 この式のすごいところは、g(x,y)=0という陰関数が与えられれば、それを各変数について陽の形にすることができなくても、変形して代入した時と同じ微分の条件が作れてしまうところです。 ラグランジュ未定乗数法は、未定乗数を使うことで、機械的に上記の計算ができるようになっています。また、未定乗数λにも意味があり、例えば今回の問題でλの値を求めると、それはVを変化させた時に、どれくらいSが変化するかを表しています。

Posted by ぼんど at 02/01/2010 05:39:38 PM

わかり易いご説明ありがとうございます。微分の公式を繙きつつ数式の変形を追って、感動致しました。 なるほど、(f-λg)'=0のλってのはdf/dgのことだったんですね。だからf+λgじゃなくてf-λgなんですか。 2つの陰関数の関係式が、1つの変数だけ陽でそれ以外は陰のままで表すことができるのがすごいですね。変数を1つだけ順に表に出してくることができる感じだと思いました。 dS/dx=0とdS/dz=0からのVの消去についても、理解しました。例えばそれらの2つの微分方程式の両方がVの2次関数だったりするとV=w(x,z)の形にするのは容易ではないのではないか、という疑問だったのですが、Vの式の%pi*x^2*yの部分を%pi*x^2*y^2や%pi*x^2*sqrt(y)にしたり、yの項を増やしたりしてMaximaで試した所、いずれも容易にV=w(x,z)の形にできました。y=f(x,z)と書けた上で本当に上の連立微分方程式からVが消去困難だと、きっとf(x,z)が陽関数とは言えない形をしているのだろうと思いました。

Posted by ynomura at 02/03/2010 09:32:31 PM

ラグランジュの未定乗数法

Lagrange multiplier EMANの物理学・解析力学・ラグランジュの未定乗数法 ラグランジュの未定乗数法 (Weblog on mebius...[link]

Posted by at 05/23/2011 08:06:45 PM

Feb 15, 2009

統計学復習メモ9: 主成分分析の計算方法

主成分分析の話になると大抵、固有ベクトルか特異値分解のどちらかが出てくる。主成分の単位ベクトルが固有ベクトルとして直ちに求められ、主成分の大きさが特異値分解を使って直ちに求められるからだ。

それにしても、行列には固有値や固有ベクトルがあるというのも、統計には分散や共分散があるというのも、理系の大卒なら誰でも知ってることだろうが、「共分散行列の固有ベクトルが主成分の単位ベクトルになる」と言われると、全く関係ない3つのものが突然シンプルに組み合わせられた感じで、衝撃的である。まるで
e^{-i¥pi}=-1
のような神秘的なものを感じるのは筆者だけだろうか。

そのこと自体は、至る所に証明というか説明があり、その方法も幾通りもある。基本的には、あるデータのど真ん中を通る直線として、最小2乗法のように2乗誤差を最小にする直線を求める代わりに、その直線方向の分散が最大になる(その直線上に射影した場合の分散が最大になる)直線を求める、という考え方である。そう言われればそういう考え方もありかな?とも思ったりするし、やっぱり最小2乗法で求める方がいいんじゃないの?とも思ったりするが、結局、主成分は分散最大で求めても最小2乗法で求めても同じ結果になるらしい。

とりあえず、固有ベクトル計算でやってみる。
m次元のデータn個を

(iは標本番号、pはパラメーター番号の意味)とし、共分散行列の固有値を

とし、対応する固有ベクトルを

とすると、W1〜Wmが第1主成分〜第m主成分、λ1〜λmがそれぞれの寄与率となる。
X,Wは行列である。プログラミング用語の「待ち行列」とかのqueueではなく、線形代数のmatrixである。フフン、行列くらい知っとるわい、という顔をしつつ、厄介だなあ、と思わざるを得ないが、主成分分析するのに行列は避けて通れない。
前回使用したデータについて、Maximaに計算させてみる。

Maximaへの入力

load(descriptive);
load(lapack);
fpprintprec:6;
X: read_matrix("data2");
[L,W,dummy]: dgeev(cov(X),true,false);
L;
W;

Maximaの出力(上のをbatch()で読み込んだ)
(%i7) L
(%o7) [.330129, .0531125, .00529463]
(%i8) W
[ .714134 .697955 - .0535787 ]
[ ]
(%o8) [ 0.4882 - .551443 - 0.67644 ]
[ ]
[ .501671 - .456912 .734546 ]

固有値(寄与率)は絶対値が一番大きいのが1つ目、二番目に大きいのが2つ目なので、第一主成分の単位ベクトルは(0.714, 0.488, 0.502)、第二主成分の単位ベクトルは(0.698, -0.551, -0.457)である。3つ目の固有値は小さいので、3つ目の固有ベクトルは無視することにする。
今回はたまたま大きい順に並んでいるが、固有値は絶対値の大きさ順に並ぶとは限らないので注意が必要である。
もしdgeev()が使えない場合は、eigens_by_jacobi()を使う。
(%i9) [L,W]:eigens_by_jacobi(cov(X)); L; W;
(%o10) [.330129, .0531125, .00529463]
[ .714134 - .697955 - .0535787 ]
[ ]
(%o11) [ 0.4882 .551443 - 0.67644 ]
[ ]
[ .501671 .456912 .734546 ]

2つ目の固有ベクトルの符号が逆だが、次に求める主成分の大きさの符号が逆になるだけなので、どちらでも良い。

次に、第一主成分、第二主成分を軸に取った時にXの各標本がどういう値になるか(主成分の大きさ)を求める。主成分の大きさYは
Y=XW
なので、

W: submatrix(W,3); /* 弟三主成分を削除 */
Y: X.W;

で求まる。これによって得られるYをグラフにすると、次のようになる。

前回のグラフと見比べると、視点が反対(前回のグラフ作成の時はengens_by_jacobi()を使ったため)だが、例えば楕円外の点の分布は大体似ていることがわかる。
これによって、今回使った3次元のデータが、情報の欠落を最小限にして、2次元で分析できる。

特異値分解による方法でもやってみる。
今度はXの中の標本は横に並べる。縦方向の1列が1つの標本である。

このXが特異値分解によって
X=U¥Sigma V^T
(Tは転置)と3つの行列(U=左特異行列、Σ=特異値行列、V=右特異行列)の積になる時、Uが主成分の単位ベクトルを縦に並べたもの、Σの対角成分(特異値)が[n×寄与率]の平方根となり、標本の主成分の大きさYは
Y=U^TX=¥Sigma V^T
となる。
Maximaで特異値分解するには、LAPACKのdgesvd()を使うと一発である。(LAPACKが使えないと、かなり厳しそうである)

Maximaへの入力

load(lapack);
fpprintprec:6;
X: transpose(read_matrix("data2"));
[s,U,V]: dgesvd(X,true,true)$
U;
s;

Maximaの出力
(%i7) U
[ - .713141 .698967 - 0.053619 ]
[ ]
(%o7) [ - .488765 - 0.55059 - .676727 ]
[ ]
[ - .502532 - .456395 .734279 ]
(%i8) s
(%o8) [5.77344, 2.30564, .727852]

(Vは100x100の行列なので省略)
行列Uに、先程と大体同じ主成分の単位ベクトルがあることと、3つ目の特異値が小さく、寄与率が低いことがわかる。

Maximaへの入力(続き)

S: zeromatrix(2,100); /* 行列Σ(3x100)の上2行分 */
S[1,1]: s[1]; /* 対角成分1 */
S[2,2]: s[2]; /* 対角成分2 */
Y: S.V;

これによって得られるYをグラフにすると、次のようになる。

さらに反対向きになったが、そういうものである。

さて、固有ベクトル計算による方法と特異値分解による方法をどう使い分けるかだが、一般に、特異値分解を使う方法の方が計算誤差が少ない、と言われる。固有ベクトル計算による方法だと、共分散行列にする時に一気に情報量が落ちるので、感覚的には納得できるが、使用するモジュールやライブラリの精度による所もあるので、一概には言えないと思う(例えばすごい桁数の共分散行列から固有ベクトルを求めれば精度は上がるはず)。

また、ライブラリの特性の問題かも知れないが、eigens_by_jacobi()やLAPACKのdgeev()では固有ベクトルが固有値の絶対値の順に並ぶとは限らないので、運が悪いと固有値を見て固有ベクトルを並び替える必要が生じる。それに対し、LAPACKのdgesvd()は特異値が必ず大きい方から並ぶので、LAPACKを使うなら特異値分解による方法の方が楽かも知れない。

See more ...

Posted at 19:15 in 数学 | WriteBacks (2)
WriteBacks

大変参考になりました

Posted by 岩淵 at 12/22/2010 05:55:32 AM

ありがとうございます。

Posted by ynomura at 12/22/2010 10:51:28 PM

Feb 14, 2009

EclipseにTomcatプラグイン導入

仕事で科学技術に恩返しできないならインターネットで、と思ってサンプルプログラムをここに書いてきたが、1年前に公開した「3択首都当てクイズ」のソースコードは掲載しなかった。なぜかというと、大体動くようになったが最後、そのコードを触るのが嫌になったからである。なぜ嫌になったかというと、Javaのサーブレットの開発が大変面倒だったからだ。

この24時間運転のwebサーバーは玄箱HGという小型の機械で動いている。いわゆるパソコンではなく、消費電力は低いが計算能力も低い。サーブレットのプログラムを書き換えると、それをこの機械で動かすためにTomcatを再起動するのに2分以上かかる。加えて、思い通りに動かないと、解析するにはTomcatのログを見るしか無い。ちょっといじって動かなくなると、調べるのも嫌になるのである。

そんな苦い思い出を、パソコンにインストールしたEclipseのTomcatブラグインが吹き飛ばした。どうせServlet APIの補完が効くのがいいくらいでしょ、と期待せずに一応インストールしてみたのだが、とんでもなく便利である。何が便利かというと、servletをデバッグモードで動かせることと、servletのソースコードを書き換えると勝手にコンパイルされてTomcatに組み込まれて動作開始することである。

Javaのサーブレットのソースコード上でブレークポイントを設定してTomcatを起動しておくと、クライアントがEclipseのデバッグモードで起動したJavaアプレットであろうが普通のwebブラウザであろうが、クライアントからの要求が来ると、サーブレットの動作がブレークポイントで一時停止する。デバッグモードのTomcatはバックグラウンドで動くので、Tomcatを起動した後も何の支障もなくEclipse上で他の作業ができる。クライアントのJavaアプレットにもブレークポイントを設定しておくと、サーブレットとアプレットとを交互に一時停止させることもできる。両方の動作を時間順に追跡できるので、大変わかりやすい。

さらに、一時停止したサーブレットのコードを、その場で書き換えて、そのメソッドの先頭からやり直させることができる。既に到着したクライアントの要求はそのままに、書き忘れた処理を追加して続行させたり、試しに処理をコメントアウトして続行させたりできるのである。ある1つの要求に対して、処理を書き換えながら色々な処理を試せるのが便利だ。
サーブレットのコードを変えたつもりなのに動作が思い通りにならない時、きちんとコンパイルし直したか、正しくTomcatに組み込んで再起動したかをいちいち気にする必要も無くなった。ブレークしてその場でコードを追加してそこを通る様子を見ればいいのである。

という訳で、EclipseのTomcatプラグインを使って見るとボロボロだった3択首都当てクイズのコードを、デバッグして整形したので、その準備過程の一部を書き残す。

●Tomcatプラグインのインストール
インストール方法は色々な所に書かれているので、省略する。注意点として、一次配布元のSysdeoのサイトは長期間アクセス不能になっているので、Eclipse Wikiのベージからダウンロードするのが良い。
今回はV3.2.1をEclipse V3.4にインストールした。

●プラグインの設定方法
EclipseのPreferencesの画面で、Tomcatの設定を埋める。
といっても、今回明示的に設定したのは、以下の項目だけ。
・Tomcatバージョン
・Tomcatホーム
・JVM設定のJREのバージョン
・Tomcatアプリケーションマネージャーのユーザー名とパスワード

●Tomcatプロジェクトの作成
取っ掛かりの部分はeclipse@Wkiのページに書いてある。
JavaのソースはWEB-INF/srcに、web.xmlはWEB-INFに、JSPのソースはプロジェクトの直下に配置する。

●web.xmlの作成
Tomcatのサンプルをコピーしてテキストエディターで書き換えても良いが、今回は、Eclipseで.xmlを開くと自動的に起動する、Eclipse組み込みのXMLエディターを試してみた。

1. EclipseでXMLファイルを新規作成(wizardはXML/XMLを使用)
2. 最初から入ってる1行("?=? xml")を消して空にする
 (これが残っているとサーブレットが起動できなかった)
3. 右クリックでAdd child→New elementを選択して"web-app"エレメント追加
4. 右クリックでEdit namespacesを選択
 Edit Schema Informationの画面でAdd→Specify New Namesp...を選択
 Prefix: (空)
 Namespace Name: http://java.sun.com/xml/ns/j2ee
 Location Hint: http://java.sun.com/xml/ns/j2ee/web-app_2_4.xsd
 を指定して追加
 これで定型部分が埋まってweb.xml用の補完が効くようになる
5. web-appを右クリック→Add attribute→versionで"version"属性追加
6. web-appを右クリック→Add Child→(略)→servletで"servlet"エレメント追加
 "servlet"エレメント以下に追加された項目を埋める
7. web-appを右クリック→Add Child→(略)→servlet-mappingで"servlet-mapping"エレメント追加
 "servlet-mapping"エレメント以下に追加された項目を埋める

●WARファイルのエクスポート
1. 作成したTomcatプロジェクトのフォルダアイコンを右クリック→Preferences選択
 →"Tomcat"の「WARエクスポート設定」タブにて.warのファイル名を設定
2. Tomcatプロジェクトのフォルダアイコンを右クリック→「Tomcatプロジェクト」→
 「プロジェクト設定に従いWARファイルを作成」を選択
3. できた.warをサーバーのTomcatのwebapps/の直下にコピー

●結果
・サーバー(Tomcat)側コード
TrichoiceServlet.java --- サーブレット本体
web.xml
QuizReader.java --- JSP用
listquiz.jsp
※ソースコード中のURL等は当方の開発用PCでのものです。

・クライアント側コード
Start.java --- メイン
MainPanel.java --- 3択画面
ServletAgent.java --- HTTP通信部分

・このサイトの起動ページ
3択アプレット
listquiz.jsp

・過去の関連エントリー
[1] [2] [3] [4]

See more ...

Posted at 16:44 in Java | WriteBacks (0)
WriteBacks

Feb 08, 2009

Maximaでの固有ベクトル計算について

筆者のMaxima 5.17+sbclの環境だと、標準のeigenパッケージを使って固有ベクトルを求めようとすると、エラーが出て求まらないことが多い。例えば、

load(eigen);
m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
eivects(m);

とやると、異様に時間がかかる上、
algsys failure: the eigenvector(s) for the1th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the2th eigenvalue will be missing.
algsys failure: the eigenvector(s) for the3th eigenvalue will be missing.

というメッセージに続いて、虚数混じりのおぞましい値が固有値として出てきて、固有ベクトルは出力されない。そもそも固有値が間違っている。

固有ベクトルを英語でeigenvectorと言うことをたまたま覚えてたので、Maximaのヘルプでeigenvectorを検索して、出てきたのがeigenパッケージのeigenvectors()/eivects()だったので、それを使ったのだが、それが間違いだったようだ。Maximaのヘルプのeigenvectorsの項にも、固有値が正しく求められないことがあると、はっきり書いてある。対処方法も書いてあるが、とても面倒である。

インターネットで調べてみると、Maximaのメーリングリストのアーカイブに辿り着いた。どうやら2007年4月くらいからこの問題は原因も含めて知られているが、まだ解決していないらしい。

数式を含まない、純粋な数値行列なら、他の方法があることがわかった。
(1) eigens_by_jacobi()
標準関数として存在する。対称行列に限るという制約があるが、例えば共分散行列の固有ベクトルを求めるなら何の問題も無い。先に言ってくれという感じである。

(%i1) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i2) eigens_by_jacobi(m);
(%o2) [[- .7932213800442011, 1.180843735860609, .5243776441835921],
[ .6946378770918785 .5975383492309946 .4005323219261915 ]
[ ]
[ - .2001963037214474 .6953728293011738 - .6902014693159894 ]]
[ ]
[ - .6909411405362675 .3992549930407618 .6026572747810012 ]

eivects()と同様、固有値の行列と固有ベクトルの行列が一緒に返される。

(2) LAPACKパッケージのdgeev()
昔からある最強の線形演算ライプラリであるLAPACKのLisp版がMaximaに同梱されている。最初に使用する時にコンパイルがなされるが、これが環境依存であり、Lispの実行環境によってはコンパイルに失敗してしまうらしい。無事にコンパイルに成功すると、dgeev()で最強の計算結果が得られる。

(%i1) load(lapack);
(%i2) m: matrix([0.123,0.456,0.789],[0.456,0.789,0],[0.789,0,0]);
(%i3) dgeev(m,true,true);
(%o3) [[- .7932213800442016, 1.180843735860611, .5243776441835922],
[ .6946378770918787 - .5975383492309948 .4005323219261915 ]
[ ]
[ - .2001963037214475 - .6953728293011737 - .6902014693159894 ],
[ ]
[ - .6909411405362672 - .3992549930407617 .6026572747810011 ]

この例は対称行列だが、もちろん非対称の行列でも問題ない。

筆者はgclとsbclで試したが、sbclではコンパイルが終了しなかったが、暴走したMaximaを強制終了すると、次からはコンパイル無しでロードできるようになった。(テスト段階で暴走していた?)

なお、root権限でインストールしたMaximaのLAPACKパッケージをコンパイルするには、root権限が必要である。(UNIX系のOSなら1回rootでMaximaを起動してload(lapack);とやれば良い)

See more ...

Posted at 23:19 in 数学 | WriteBacks (2)
WriteBacks

eigen_by_jacobi関数、役に立ちました。 ありがとうございます。 一つ質問があるんですが、数式を含んだ状態で固有値を求めた際、虚数を含まないようにする命令はありますか? A:matrix([-z,-1,0],[-1,-z/2,-0.25],[0,-0.25,-z/2]); B:eigenvalues(A); として固有値を求めると、対称行列なのに虚数が表われて、よく分からない数式になってしまいます。

Posted by yuji at 12/07/2010 12:51:00 PM

eigenvalues()はeigenパッケージにあり、まだeigenvectors()と共にバグがありますので、不正な答えが出力されることがあります。 また、3x3の大きさの行列の固有値は3次方程式の解であり、3次方程式の解の公式は(実数解でも一時的に)虚数が含まれますので、行列に変数が残っていると、実数解しかあり得ないとわかってても、虚数抜きで表現するのは困難だと思われます。(ちなみにMathematicaだと、こういう場合は「カルダノの公式」が適用されて、やはり実数解が虚数を使って表現されるようです) 固有値は固有方程式の解であることを使って、Maximaで solve(determinant(A-diagmatrix(3,x))=0,x); とすると、虚数は入りますが、正しい答えが得られます。 もし、zに何か実数が代入されてから計算されても良いのであれば、 A(z):=matrix([-z,-1,0],[-1,-z/2,-0.25],[0,-0.25,-z/2]); B(z):=eigens_by_jacobi(A(z)); としてから B(3) などとすると、実数表示の解が得られます。(ここでも残念ながら、eigenvalues()を使うと虚数を含む解の公式表示になってしまうようです)

Posted by ynomura at 12/07/2010 10:59:57 PM

Feb 06, 2009

統計学復習メモ8: 回帰直線と主成分分析

高校生の時に、必要も無いのに関数電卓を買った覚えがある。物理の先生に、理系ならずっと使うから早目に買っといて損は無いとか言われたからだと思うが、結局大学で物理実験をやるまでは遊びでしか使わなかった。それでも、元々電卓が好きだった私は、毎日毎日、関数電卓のプラスチックのボタンをカチャカチャと押して遊んでいた。
その電卓には、「一次回帰直線」の係数を求める機能があった。学校では習わなかったが、説明書にy=a+bxのaとbを求める機能と書いてあったので、意味は分かった。来る日も来る日も自室に逃げ込んで電卓と戯れてる内に、「一次回帰直線」がとても身近なものになった。

そんなに一次回帰直線に慣れ親しんだから、主成分分析には敏感に反応した。というのも、主成分分析は繰り返し回帰直線を求めることに似ていると思うからである。

一次回帰直線とは、パラメトリックな2変量の関係を直線で近似するものなので、2次元のデータを1次元で近似することに相当する。2変量をx,yとすると、標本(x1,y1)〜(xn,yn)の関係を直線y=a+bx(a,bは定数)で近似するので、その直線との誤差を無視して(x,y)=(0,a)+(1,b)tと置くと、2変量の標本(xi,yi)が1変量tiで表せるのである。
主成分分析の第一主成分だけを使うことを考えると、多変量の関係を1次元で近似することになるから、k次元の標本(x1i,x2i,…,xki)を(C1,C2,…,Ck)+(p1,p2,…,pk)tiに近似すると、k個の変量を持つ標本(x1i,x2i,…,xki)が1変量tiだけで表せる。(C1,C2,…,Ck)は定数ベクトルで、(p1,p2,…,pk)が第一主成分の単位ベクトルである。第二主成分の単位ベクトルを(q1,q2,…,qk)としてこれも使うとすると、(x1i,x2i,…,xki)≒(C1,C2,…,Ck)+(p1,p2,…,pk)si+(q1,q2,…,qk)tiという感じに、k変量の標本を2変量で近似することになる。

一次回帰直線の最小2乗法による求め方をおさらいする。標本(xi,yi)と直線y=a+bxとの誤差はyi-a-bxiであるが、この誤差の2乗和を最小にするようにa,bを決めたい。誤差の2乗和をEとすると、E=Σ(yi-a-bxi)2で、このEはaの2次関数であり、bの2次関数でもあるから、Eをaで偏微分したものが0になるaと、Eをbで偏微分したものが0になるbを求めれば良い。
\frac{\partial E}{\partial a}=0, \frac{\partial E}{\partial b}=0
これを解くと、
b\frac{\sum xy - \frac{1}{n}\sum x\sum y}{\sum x^2-\frac{1}{n}(\sum x)^2},a=\bar{y} - b\bar{x}
が得られる。bの分子はxとyの共分散、分母はxの分散なので、xの分散をVar(X)、xとyの共分散をCov(X,Y)と表すと、b=Cov(X,Y)/Var(X)とも書ける。実はシンプルである。

念のため、試しにやってみて確認する。例によってMaximaを使う。サンプルデータをdata1とし、

data1: read_nested_list("data1");
datam: read_matrix("data1");
load(descriptive);
b: cov(datam)[1][2] / var(datam)[1];
a: mean(datam)[2] - b * mean(datam)[1];

とすると、b=-.8338664607943058、a=5.848350705139143と求まる。そこで
plot2d([[discrete, data1], a+b*x],[x,0,2],[style,[points],[lines]]);

とすると、次のグラフが作られる。

主成分分析についても、詳細を省いて、とりあえずサンプルデータ(data2)を使って図示してみる。
下のように、標本が原点を中心に楕円盤状に散らばっているとする。(どれかの画像をクリックするとgifアニメが開きます)



少しずつ回転させてみたが、わかりづらいので、大体の分布の形を枠線で表示する。


第一主成分の単位ベクトルは最も誤差の和を小さくする直線の方向、すなわち楕円の長径方向の長さ1のベクトルになる。次の図の青線が第一主成分(をスカラー倍したもの)である。


標本から第一主成分を抜く(第一主成分を0にする)と、下の図のようになる。(回転方向は変更)


つまり、第一主成分と垂直な平面への射影である。元が楕円盤状だったのでこれも楕円状であり、この楕円を最も近似した直線、すなわち長径方向が第二主成分の向きになる。(下図緑線)


3次元なので、第二主成分も抜くと、直線状になる。その直線は第一主成分、第二主成分と直交し、その直線上に第三主成分がある。

空間が何次元でも同様である。最も誤差が小さくなる近似直線の方向の大きさ1のベクトルが第一主成分の単位ベクトルであり、その方向のベクトル成分を抜いた(超平面に射影した)ものの近似直線の方向の大きさ1のベクトルが第二主成分の単位ベクトルであり、その方向のベクトル成分も抜いて…を繰り返すと、全ての主成分が出てくる。

そんな反復計算をせずに主成分分析する方法は、次回まとめる。

See more ...

Posted at 22:30 in 数学 | WriteBacks (0)
WriteBacks