« MIT Python 教室 第2回 (演算子、命令、条件分岐、繰り返し) | トップページ | Python Copy Special Exercise in Japanese »

2010-03-21

MIT Python 教室 第21回 (シミュレーション結果の確認、曲線のあてはめ、線形回帰、PyLab)

Python と計算機科学の入門講義の21回目.

前回は円周率をモンテカルロ法により求める方法を説明していました. 今回は, 大きく分けると前半と後半に分かれています. 前半は, その求められた円周率の値の正確さを議論し, 推定した値の正しさは, 推定に使ったモデルに依存することを説明しています. そこから発展して, 後半では 測定したデータに対するモデル (関数) の当てはめ (最小二乗法) について説明しています. 最小二乗法について, 計算方法の導出などの説明はありませんが, Python のライブラリである Pylab に含まれる最小二乗法をするための関数を説明しています. データに対する関数の当てはめについては, 次回の講義に続きます.

まず, 黒板に 3.14161124 という 数字が書かれているところから講義がはじまります. これは, モンテカルロ法により 10億個の点で推定された円周率とのことです. 100億個の点で推定した値もほぼ同じになりました. では, どれだけの点が推定に必要なのでしょうか? そして, どのくらい正確な値が得られるのでしょうか? この2つの問いは非常に関係しています.

例えば, MIT にいる中国人の学生数を 推定するとき, 学生を数百人調査すれば通常は十分ですが, 均一でない調査になった場合には 正確でない推定になってしまいます. 本当に正確な推定を行うには, 母集団全体を調査する必要があります.

さらに, 一般には推定した値が正しいかどうかを 確認する方法はなく, その値がどのくらい正しそうなのかを 知ることが出来るだけです. 円周率の推定の場合, 何回も繰り返して求めた推定値がどれもほぼ同じ値に 収束すると, ほぼ正しそうだと考えられます.

しかし, 同じ値に収束するからといって正しい値とは 限りません. もし, 円周率を求めるときに使った計算式の 係数を 4 から 2 に変えても, 値は 1.57 付近に常に収束しますが, その値はあきらかに円周率ではありません.

この例のように, 値を推定するための モデル自体に誤りがあれば, 収束している値が推定値として間違っていることは ありえるのです. ある値に収束しているという統計的な妥当性は 必ずしも推定値の正しさを保証しません.

モデルが間違っていることをどうやって 見つければよいか, という問題が出てきますが, 実際の例に当てはめて確認することが必要になります. つまり, 推定値の正しさを確認するには, 実験をして求めたデータとの整合性を見る必要があります. ということで後半に入っていきます.

データに対して, そのデータを 説明するモデルを考え, その モデルから結果/予測を求める, と いうのが, 科学や工学の仕組みであると説明しています.

ここでは, 物理の古典的な問題として, ばね定数を求める問題を考えます. あるばねは伸ばしやすく, またあるばねは伸ばすのに力が必要になる, という違いは, それらのばねの ばね定数に違いがあるためです. フックの法則によると, ばねに かかる力 F とばねの伸び x に対して,

F = -kx
という比例関係が成り立ち, その中の k がばね定数にあたります. ばねの伸びに対して, 力が負の値に なっているのは, 伸びの方向と ばねが戻ろうとする力の方向が逆に なっているためです.

ではどうやって, ばね定数を求めればよいのでしょうか? ここでは, ばねの伸びと掛かる力のデータから ばね定数を求めます. そのプログラムでは, "距離:力" という書式のデータを 読み込んでいますが, 読み込む方法の詳細は説明されていません. 掛かる力に対するばねの伸びのプロットを 表示しています. フックの法則によれば, 点は一直線に並ぶはずですが, プロットでは直線に沿う形で 点が広がっています. 測定データが直線にならないのは, 測定誤差などが入っているためで, このようなノイズが含まれるのは 避けられません. その点で, 実際のデータは 理論にぴったりと合うことはないといえます.

そこで, 測定データに直線を当てはめ, その直線の傾きをばね定数として求めます. なぜ直線なのか, 2次曲線や円ではないのか, という疑問に対しては,

  • 1. モデル (フックの法則) が直線を薦めている,
  • 2. 測定データが直線のように並んでいる,
ことによります.

直線を当てはめるには, 最適な直線を見つける ための目的関数が必要になります. ここでは, 最小二乗法に基づく目的関数を 使います. 観測値と予測値の二乗誤差の和を 最小化するのが目的になります. ここで誤差を二乗するのは, 正負の誤差が 互いに打ち消しあうことを避けるためです.

直線を見つける方法として, 単純に目的関数の値を見ながら 傾きを変化させたり, 以前の講義で 見たような Newton-Raphson 法を使うことも 考えられますが, 最小二乗法の場合は 閉形式で解を求めることができます.

ここでは, pylab に備え付けの関数 polyfit() を使います. polyfit() は, 点の集合が与えられたとき, それらの点に対する最小二乗近似となる 多項式関数を見つけます. ばねの伸び distances, ばねに掛かる力 forces に対して,

k, b = polyfit(distances, forces, 1)
は最小二乗法により求めた直線 y = k * x + b の係数 k, b を返します. polyfit() の最後の引数 1 は, 当てはめる多項式が1次式 (直線) で あることを示しています.

この手順により見つかった k, b により直線をグラフに表示しています. 観測データ中の点は直線上にはあまりなく, 大きく離れた外れ値もありますが, 多くの点が直線の近くにあります. これは二乗誤差を最小化したことによります.

さらに, 別の実データを使った例を考えます. 時間に対する速度の変化をグラフに表示して, 直線を当てはめます. そうすると, あまりこのデータには直線は当てはまって いないようです. データ点の並びを見ると, 放物線のようなので, 2次曲線を当てはめます. コマンドを使った2次曲線の当てはめは,

a, b, c = polyfit(times, speeds, 2)
とします. 直線より, 2次曲線の方が当てはまっているように 見えます.

ここまで行ってきたことは, 線形回帰と呼ばれます. 従属変数 y と 独立変数 x との 関係が y= ak xk +...+ a2 x2+ a1x +b となっていて, 係数 ak,..., a1 が1次である式を, 線形回帰は対象とします. x は多項式になっても構いません. (当てはめを行うときは, x の項は具体的な 値が代入されるので, 当てはめの困難さに影響しない) 今回は, 変数は x のみになっています.

黄色で図示された2次曲線よりも, 赤色で図示された直線の方が当てはまりが 良くないように見えるというのは, 主観的な見方に過ぎません. 当てはまりの良さを客観的に議論するために, 形式的な方法を導入します.

決定係数 R2 R2=1-(EE/DV) と定義されます. ただし, EE は推定における誤差, DV はデータ中の分散をそれぞれ表し, EE=\sumi(yi-fi), DV=\sumi(yi-yave) 2 により計算されます. yi が観測された従属変数の値, fi が予測された従属変数の値, yave が観測された従属変数の平均値 とします. R2 = 0.9 ならば, データの 90% は推定されたモデルで 説明され, 10% は誤りであることを 示しています. R2 が 1 より大幅に小さいときには, 当てはめが十分でないと考えられます. 逆に当てはめが完全なとき, つまり, すべての観測点が予測点と一致するなら, R2=1 になります.

時間に対する速度の変化について, 直線の当てはめに対する決定係数は R2 =0.896, 2次曲線の当てはめに対する決定係数は R2 =0.973 となりました.

2次曲線の当てはめの方が 直線の当てはめよりも良い決定係数の値が得られました. これは一般的にもそうで, 決定係数の値は 高次式を当てはめたときほど1に近づきます. その理由は, 高次式 (例えば2次式) は 高次元項の係数を 0 (a2=0) とすることで, 低次式 (1次式) で表すことが できるためです.

7つの点が与えられたとき, 直線と5次関数を当てはめたときの 決定係数の値はそれぞれ 0.978, 0.99 になりました. しかし, 良い決定係数の値が 良い当てはめになっているわけではありません.

フックの法則に対するデータに 5次式を当てはめた場合を示し, ばねの伸びが十分大きいとき, y の値が x の増加の比よりも大きく 予測されることを示しています. これは, 当てはめた関数が データに適合しすぎているためです. より複雑な関数になるほど, 過適合になる可能性が増えていきます. この問題についての説明は次回に続きます.

« MIT Python 教室 第2回 (演算子、命令、条件分岐、繰り返し) | トップページ | Python Copy Special Exercise in Japanese »

Python lecture notes」カテゴリの記事

コメント

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.f.cocolog-nifty.com/t/trackback/1174621/33830727

この記事へのトラックバック一覧です: MIT Python 教室 第21回 (シミュレーション結果の確認、曲線のあてはめ、線形回帰、PyLab):

« MIT Python 教室 第2回 (演算子、命令、条件分岐、繰り返し) | トップページ | Python Copy Special Exercise in Japanese »