ラグランジュ方程式のイメージ解釈
ラグランジュ方程式とは、解析力学を初めとする物理学の基礎を成す方程式です。
∂L( q(t), q'(t), t )/∂q(t) - d/dt{∂L( q(t), q'(t), t ) /∂q'(t) } = 0
-- wikipedia:オイラー＝ラグランジュ方程式
ちょうど関数の最小値を求めたいとき「微分＝０」を調べるのと同じように、
汎関数（関数の関数）の最小値を求めたいとき、このラグランジュ方程式を調べます。
ただ、「微分＝０」には“谷底の傾きは平らになる”という明確なイメージが描けるのに比べて、
ラグランジュ方程式のイメージを思い描くのはかなり難しい。
数学が相当得意な人であっても、記号の字面を追うのが精一杯、というのが実情でしょう。
ラグランジュの『解析力学』という本には、イメージを助ける図が一切ありません。
「緒言」でラグランジュ自身が述べているように，この本には「図がまったく見出され」ず，必要とされるのは「代数的な操作のみ」である．
-- <翻訳> J・L・ラグランジュ『解析力学』(抄)
-- https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/137415/1/phs_5_127.pdf
そもそも解析力学とは、
「代数的な操作こそが正当であり、幾何学的な図は本来排除されるべき補助手段に過ぎない」
という思想の産物だったのです。
この思想は現代でもなお有効で、高度抽象的な対象を前にして、私たちは
イメージ化をスッパリ諦めるか、
邪道であっても図を思い描くか、
の２択を迫られているわけです。
ラグランジュ方程式が難しいと感じる本当の理由もここにあって、
最初からイメージ化を諦めてかかった方がむしろ受け容れやすい、という側面があります。
にもかかわらず、ここではあえて邪道なイメージ化の方法を試みます。
ラグランジュさんには申し訳ないのですが、私たち凡人はイメージ化して初めて腑に落ちたと感じるからです。
ラグランジュ方程式の骨格を読み解くため、いったん“微分”という要素を棚上げし、
まず“関数の関数”という要素だけに着目しましょう。
いま、
L(t) = t^2 - 2 t
という、t の関数があったとしましょう。
ここで L の最小値を求めよと問われれば、関数 L(t) を t で微分するのが常套手段なのですが、
ここではあえて次のような見方をとります。
T(t) = t^2
U(t) = 2 t
と置いて、
L( T(t), U(t) ) = T(t) - U(t)
という“T と U の関数L”を作ります。
このとき L の最小値を求めるには、
T と U がちょうどバランスをとって等しくなる点、T' = U' を探せばよい、
ということになるでしょう。
T'(t) = 2 t
L'(t) = 2
2 t = 2 を解いて、t = 1 のとき L は最小となる。
こんな手間をかけて何が嬉しいのかというと、最小値についてのイメージが変わります。
関数 L(t) を直接 t で微分するとき、最小値とは“傾きが０となる点を探す”ということでした。
一方、L(T,U) = T - U を間接的に t で微分したとき、最小値とは“T' と U'がバランスをとって打ち消し合ったとき”という意味に変わります。
“傾きが０となる点”という見方をするなら、最小値は“点”として得られました。
一方、“バランスをとって打ち消し合う”という見方をするなら、
T'=U' を満たす点を列挙することで、最小値を“線”として得る道が拓けます。
ここに変分法 ～ 最小となる関数を求める方法 ～ のカラクリがあります。
解析力学における「ダランベールの原理」も、結局はこの“バランスをとって打ち消し合う”ことを意味しています。
運動する物体は、局所局所において、静的な力と、動的な加速（に由来する力）が打ち消し合っている。
なので、そうしたバランスを保つ点を線で結べば、運動の軌跡が自ずと浮かび上がるだろう。
これが、ラグランジュ方程式を形作るアイデアだったのです。
“バランスをとって打ち消し合う”ことを経済に例えるなら、需要と供給の均衡点のようなものです。
-- wikipedia:需要と供給 より引用
ラグランジアンとは、余計なエネルギーコストという意味で、経済で言うところのコストと高い親和性があります。
仮に「需給のバランスがとれた姿が実現すべき社会である」と考えるなら、
需給バランスを原理原則として、社会の姿を予言（計算）できることになるでしょう。
たとえば原料価格を横軸にとって、均衡点を線で結べば、原料価格に対する適正な生産量のグラフが描けます。
あるいはもし、消費者の満足度のようなものが数値化できれば、
満足度を横軸にとって均衡点を線で結べば、適正価格のグラフが作れるはずです。
※ 私は経済学には疎いのですが、どうも昨今の経済学（の一派）は、そういう事をやっているみたいです。
※ 『大学二回生のときに「解析力学」の講義で習ったラグランジュの方程式が、なぜかミクロ経済学の教科書にそのまま書かれているのである。』
※ * 現実と違っても数学モデルを信仰してしまう経済学者たち >> http://data11.web.fc2.com/jiyuuboueki5.html
需給バランスの例で言うと、オリジナルのグラフは横軸が「数量」となっています。
この１枚のグラフから得られる答は「均衡点」というただ１点のみであり、
それだけで終わっていたなら、需給バランスのアイデアはここまでありがたがられなかったことでしょう。
需給バランスが有用なのは、オリジナルのグラフの上にさらにもう１次元、原料価格とか、満足度とかいった、新たな変数が追加できるところにあります。
オリジナルの均衡“点”は、新たな変数 ～ 原料価格に対する“線”、あるいは満足度に対する“線”となります。
さて、だいぶ持って回した言い方となりましたが、先ほどの単純な例、
L(T,U) = T(t) - U(t) = t^2 - 2 t
に立ち返ると、そこにあったのは変数 t のみで、このままで終わっていては何のありがたみも感じられません。
変数 t の上に、新たな変数、位置 q と、速度 q' をオーバーラップさせて初めて、
・位置 q に対する均衡の線
・速度 q' に対する均衡の線
が描けるわけで、ラグランジュ方程式の真価が発揮されるのです。
※ 解析力学では、物体の位置を x と書かずに q と書くことが多いです。
※ x がいわゆる直交座標を表すのに対し、q は広く一般的な座標という気持ちが込められています。
改めて、先ほど省いてきた“微分”の要素に目を向けましょう。
ラグランジアンＬを、本来の姿、
L( q(t), q'(t) )
に戻します。
つまりラグランジアン L とは、位置 q と、速度 q' の関数です。
※ 本来であれば、さらに時刻 t も含めるべきなのですが、簡単のため省略します。
冒頭に掲げたラグランジュ方程式の骨格のみを抜き出すと、こうなります。
∂L/∂q = d/dt{ ∂L/∂q' }
ラグランジアン L を「エネルギーコスト」と読み替えると、この式は、
「エネルギーコスト L に対する、位置 q の寄与と、速度 q' の寄与が均衡している」と読むことができます。
あるいは、「静的な力と動的な力が均衡している」と読むこともできるでしょう。
（そして、位置と速度が均衡しているとき、エネルギーコストは局所的に最小となる、というのがこの式の使い方です。）
ただ、式を見ると、速度 q' のある右辺には d/dt という操作が施されています。
この d/dt を、どう解釈すべきでしょうか。
最も邪道にして簡単な解釈は“約分”です。
速度 q' という記号を、ライプニッツ流の書き方 dq/dt に直してみましょう。
（どちらの書き方でも、t で微分する、という意味は同じです。）
そして、d/dt という記号は“約分”できる、とするなら、
このように、ラグランジュ方程式とは（少なくとも強引な字面の上では）極めてアタリマエのことを主張していたのです。
“約分”よりもう少しまっとうな解釈は、「単位時間当たり」という操作の必然性です。
d/dt には「単位時間当たり」という意味があります。
仮に、時間がゆっくり流れているスローモーションの世界があったとしましょう。
スローモーションの世界では、時間はゆっくりになりますが、運動の軌跡は同じ形のままです。
時間がゆっくりになっても、位置 q の寄与（位置 q に依存する静的な力）は変わりません。
一方、時間がゆっくりになったとしたら、速度 q' の寄与（速度 q' に依存する動的な力）は小さくなるはずです。
それでも軌跡が同じだということは、速度 q' の寄与を、時間の流れる速さに合わせて調整しなければなりません。
どのように調整するかというと、時間の流れる速さで割ってやる、つまり「単位時間当たり」の量に換算します。
これが速度 q' の寄与 {∂L/∂q'} に d/dt という操作を施す理由です。
もし時間がゆっくり流れる世界が認められないのなら、時間の単位を変えて計算することを考えてみてください。
運動の軌跡は、秒単位で計算しようとも、分単位で計算しようとも一致するはずです。
にもかかわらず、単に {∂L/∂q'} を計算すると、秒速の方が分速の 60倍となり、静的な力 ∂L/∂q と釣り合いがとれなくなってしまいます。
（静的な力の方には、陽に単位時間が入っていないことに注意。）
方程式が成立するためには（つまり方程式の両辺の単位を合わせるには）、{∂L/∂q'} を時間の単位で割る必要があります。
※ この時間単位を揃える操作 d/dt と、そもそも速度とは位置の時間変化であった dq/dt とは、表記の上では“約分”のように打ち消し合います。
以上のお膳立てが整ったところで、ラグランジュ方程式のイメージ化を試みましょう。
まず横軸に時刻ｔ、縦軸に物体の位置ｑをとり、運動の軌跡を描きます。
これに高さの軸Ｌを追記し、３次元のグラフを形作ります。
Ｌは物体の位置ｑに依存するので、それをｑ×Ｌ平面上に描きます。
ｑ×Ｌ平面上に描かれた軌跡の傾きは、∂L/∂q を意味します。
ここまでは静的な力の話で、図の上では青色で描かれています。
次に、運動の軌跡を微分した速度ｑ’の線をｔ×ｑ平面上に重ねて描きます。
速度ｑ’の軸は、位置ｑの軸とスケール（単位目盛り）が異なっているため、改めて図の右側にｑ’軸を付け加えました。
つまり微分したｑ’の線は、右側のｑ’軸の方を用いたｔ×ｑ’平面上にあるものと考えます。
（図の上で、ｔ×ｑ平面とｔ×ｑ’平面は重ねて描かれています。）
Ｌは物体の速度ｑ’にも依存するので、それを描く平面が必要です。
図の上に、速度ｑ’軸に直行する形で、高さ軸 d/dt(Ｌ) を付け加えましょう。
（図の上で、右端で垂直に立っている平面がｑ’×dt/d(Ｌ)平面となります。）
この ｑ’×dt/d(Ｌ)平面上に描かれた軌跡の傾きは、d/dt(∂L/∂q') を意味します。
ここまでが動的な力の話で、図の上では赤色で描かれています。
こうして描いた３次元のグラフを、横方向から眺めてみましょう。
横軸にｑとｑ’のオーバーラップ、縦軸にＬとdt/d(Ｌ)がオーバーラップした、
２枚の青と赤の重なったグラフが見えることになります。
青いグラフの上には、Ｌに対する位置 q の寄与（位置 q に依存する静的な力）が描かれています。
赤いグラフの上には、Ｌに対する速度 q' の寄与（速度 q' に依存する動的な力）が描かれています。
このとき、Ｌを局所的に最小とするのは、２つの力の均衡点となります。
逆に、この均衡点を計算しながら辿っていくことで、もとの運動の軌跡が再現できます。
そして、均衡点を表す式とは（静的な力）＝（動的な力）という形をとり、記号で書けば
∂L/∂q = d/dt{ ∂L/∂q' }
これがラグランジュ方程式です。
* まとめ *
・微分というツールを、“傾き＝０”ではなく、“２つの勢力がバランスをとって打ち消し合った均衡”であると捉える。
・ラグランジアンＬを、「エネルギーコスト」と読み替える。
・２つの勢力を、
(1) エネルギーコストＬに対する位置 q の寄与（位置 q に依存する静的な力）
(2) エネルギーコストＬに対する速度 q' の寄与（速度 q' に依存する動的な力）
に対応させる。
・(1) と (2) は同じ土俵（座標系）の上には乗らない。
なぜなら、静的な力が単位時間を含まないのに対し、動的な力は単位時間を含むからである。
そこで、動的な力の方に d/dt という操作を加え、同じ土俵に乗せる。
・こうしてできた均衡の式が、ラグランジュ方程式だ。
∂L/∂q = d/dt{ ∂L/∂q' }
最後に、以上のようなイメージ解釈は冒頭にも述べた通り、ラグランジュさんの精神からすれば邪道なのであって、
いわゆる教科書流の「代数的な操作のみ」による証明が正道であるとされています。
しかしなぜ、代数的操作は正しく、イメージ解釈は受け入れがたいのでしょうか。
天才ラグランジュさんとて人の子であり、おそらく誰よりも明確なイメージを抱いていたはずです。
そのイメージを公開しなかった理由の１つは、逆説的ですが、
『イメージは正しく人に伝わらないから』
ではないかと思うのです。
凡才である私がラグランジアンの図を描きつつ浮かんだのは、
「この図の意図は、なかなか人には伝わらないだろうな・・・」
という思いでした。
図の読み手が悪いという意味ではなく、この図では意図が描き切れていない、表現力が足りないという意味です。
もう少し具体的に言うと、“時間の流れを揃える”といった感覚を図にすることができません。
もちろん、自然言語で表すことも難しい。
インターラクティブなアニメーション、dtをビヨーンと引き延ばすと、それに伴って {∂L/∂q'} が変化するようなオブジェクトを作れば、もう少し表現できるかもしれません。
それでも完全ではありません。
勝手な想像ですが、天才ラグランジュさんは、自分の持つイメージが当時の人たちにはまず伝わらないと悟ったのではないでしょうか。
その上で、驚くべき方法 ～ 代数的操作のみによって表現する方法を開拓したのだと思うのです。
それがあまりにもすごかったものだから、後世の人がすっかり「代数的な操作のみが真実」であると信じ込んでしまった。
確かに、代数的な記号操作はイメージよりも誤解が少ないし、コピーするのも難しくはありません。
しかし、ただコピーしただけでラグランジュさんの意図が伝わったのかというと、そこには依然、大きなギャップがあるものと思うのです。
代数的操作による説明は、以下を参照。
* ラグランジアンに意味は無い >> d:id:rikunora:20090327
平均と分散の物理モデル
以下は、前エントリー「なぜ分散は２乗の和なのか」d:id:rikunora:20190409 から出たオマケの与太話。
なぜ分散は２乗なのか、と問われたとき、「分散を最も小さくする点が平均値だから」というのが１つの答でした。
この答を物理的なモデルに当てはめると、こんな絵が描けます。
・各データと、とある１点をバネで結んだとき、１点は平均点で止まる。
このモデルで分散は、「バネに溜まったエネルギー」に対応付けられます。
分散が足し算できること（加法性）は、「エネルギー保存則」として理解できます。
ところが分散は、こんな物理モデルに当てはめることもできます。
・各データに重みを付けてバランスをとったとき、重心は平均点と一致する。
このモデルの場合、分散は「慣性モーメント」に対応付けられます。
分散が足し算できること（加法性）は、「角運動量保存則」として理解できます。
以下、確率・統計入門（小針あき宏）という本からの引用。
これだけを見ていると何でもないようだが、これは‘平均の持つ意味’についての、一つの思想を物語っている。
・・・平均値をその分布の代表とする理由は、そのまわりの分散が最小ということによって特徴づけられているわけである。
ρ(x)を質量密度と考えると、分散は慣性能率に対応していることは、説明を待たないであろう。
そして自由な剛体の運動は、加えられた力をＦ，質量をＭとするとき、
Ｆ＝Ｍｘ''
によって記述される、重心の曲線運動と、
慣性能率をＩ，重心からの力の加えられたところまでの距離をｒとすると
Ｆｒ＝Ｉθ''
によって記述される、重心回りの回転運動によって完全に決定されるが、
どうして《重心のまわり》に回転がおこるか、というと、そのまわりの慣性能率が最小だから、に他ならない。
そして重心～平均、慣性能率～分散、の対応を見るとき、平均の持つ意味について、新たな理解が生ずるだろう。
エネルギーとは、「物理的に、時間変化に対して不変な保存量」のことです。
角運動量とは、「物理的に、角度の変化に対して不変な保存量」のことです。
物理学には“ネーターの定理”と名付けられた定理があります。
標語的に言うと、「対称性あるところに保存則がある」といった主張です。
ネーターの定理によると、
・時間並進の対称性から「エネルギー保存則」が、
・回転の対称性から「角運動量保存則」が、
・空間並進の対称性から「運動量保存則」が、
それぞれ導かれます。
ならば、分散を「運動量保存則」に結びつけるような物理モデルは無いものだろうか・・・
少し考えてみたのですが、うまいモデルは思いつきませんでした。
それもそのはず、エネルギーと角運動量が２次式なのに対し、運動量保存則は１次式なので、土台あてはまらないのです。
※ 運動量保存則は、”分散”ではなく”平均”に結びつけられます。
※ 平均は、確かに平衡移動について対称、一斉にゲタを履かせても不変ですが、あまりにもアタリマエなので意識されません。
話はこれだけなのですが、以上は何も分散だけに限った話ではなく、２次式で表されるような量全般に言えるように思えます。
もともと２次式に当てはまるような量であれば何であれ、
エネルギーや角運動量に関連付いたもっともらしい物理モデルとして表される、ということです。
一般に、理解の仕方や解釈は、個々人の感性やセンス、才能といったものに委ねられると言われがちです。
しかし、この分散の物理モデルを見て思うことは、
「思考の筋道とは全くの自由奔放ではなくて、自ずと規定されている」ということです。
不自然さを感じさせない発想は、文字通り“物の理”に叶っているのです。
そうした“物の理”に叶った発想だけが、分かりやすく、人に伝えやすく、結果的に生き残るのではないか。
そんな風に思えるのです。
- 作者: 小針アキ宏
- 出版社/メーカー: 岩波書店
- 発売日: 1973/05/31
- メディア: 単行本
- 購入: 5人 クリック: 130回
- この商品を含むブログ (23件) を見る
なぜ分散は２乗の和なのか
Ｑ．なぜ分散は、単純な差（偏差の絶対値）ではなく、差の２乗を計算するのか？
Ａ．分散を最も小さくする点が平均値だから。（単純な差を最も小さくする点は中央値となる。）
“分散”というキーワードは統計学の基礎中の基礎であり、どんな教科書にも“平均”の次くらいに載っていることがらです。
しかしながら、いきなり登場する“分散”の意味が分からず、統計学の入り口で挫折する人は少なくありません。
偏差の２乗の平均、つまり、各値と平均との差の２乗の平均を分散といい、
分散の平方根の正の方を標準偏差という。
統計で、ちらばりを表すものとして、標準偏差や分散が多く用いられる。
-- 高校の教科書（啓林館）より.
教科書にはこのように書かれているのですが、これで分かった気になるでしょうか。
・なぜ、差の２乗を計算するのか？
・差そのものであってはいけないのか？
・なぜ、分散と標準偏差の２種類があるのか？
最後の疑問については「分散が２乗した結果なので、元の単位に戻すために平方根をとる必要がある」のだと説明されます。
(標準偏差) ＝ √(分散)
だとすると、そもそも値を２種類作らなければならなかった理由は２乗を計算するからであって、
最初から差そのものを採用していれば１種類で済んだはずです。
各データと平均値との差のことを「偏差」と言います。
分散とは、偏差を２乗した値の平均のことです。
一方、偏差の値そのもの（絶対値）の平均は「平均偏差」と呼ばれています。
分散（と標準偏差）に比べて、平均偏差はあまり有名ではありませんし、全ての教科書に載っているわけでもありません。
もし何の予備知識も無しに「ばらつきを数字で示せ」と言われたなら、
データの差そのものを持ってくる平均偏差の方がよほど自然に思えます。
なのになぜ、今日の統計では、ややこしい分散（と標準偏差）の方が主流なのでしょうか。
歴史を振り返ってみると、偏差の２乗が一般化する以前に、差そのものをばらつきの指標と捉えていた時期がありました。
例えばラプラスは、差そのものを誤差損失の指標として扱っていたという経緯があります。
はっきりと２乗を指標とするようになったのは、ガウスの登場以降のことなのです。
Laplaceはこのことを同じような方法で考察した。しかし彼は、つねに正として扱われる誤差自身を損失の量として選んだ。
-- ガウス誤差論より.
■ 最小２乗法のこころ
ではなぜ、ガウスは２乗を考えたのか。
その心を知るには、ガウスのお家芸である“最小２乗法”をひもとくのが早道です。
最小２乗法とは、誤差を含んだデータに対して、もっともよく当てはまる関係を見出す方法です。
例えば、とあるグループの身長と体重のデータが、こんな風だったとしましょう。
ここで「身長が高いほど体重が重い」という関係を１本の直線で示すには、どのように線を引くのがもっとも合理的でしょうか。
最小２乗法では、次のように考えます。
一本の棒と、それぞれのデータの間をバネで結んで引っ張ると、棒はほどよいところに落ち着きます。
この落ち着いた先の状態を「もっとも当てはまりが良い」と見なすのが、最小２乗法の考え方です。
なぜ２乗なのかというと、バネは引っ張れば引っ張るほど、長さに比例した力がかかるからです。
バネの長さを２倍に引き伸ばすには、
・長さが２倍
・力も２倍
合わせて２×２＝４倍のエネルギーが必要です。
これが２乗の由来です。
最小２乗法とは、つまり「バネに溜まったエネルギーが最小となる位置を探す方法」のことだったのです。
* ∫ニョロっと伸びるはエネルギー、バネは答を知っている
>> http://miku.motion.ne.jp/multi/01_MinSquare.html
いま、２次元のグラフで見た方法を、さらに単純化して１次元としてみましょう。
複数個のデータと、自由に動かせる、とある１点をバネで結んだとしたら、とある１点はどこに落ち着くか。
落ち着く先は、ちょうどデータの平均値となります。
※ ２次元の場合、当てはめた直線（回帰直線）は必ず平均値を通ります。
平均値とは、「バネに溜まったエネルギーが最小となる点」のことだったのです。
そして「バネに溜まったエネルギー」とは、各データと平均との差の２乗、つまり分散に比例しています。
かくして手元のデータから、もっとも当てはまりの良い点を探そうと考えたとき、分散の最小化というアイデアに導かれるわけです。
■ 分散の加法性
ここで１つクイズ。
『厚さが平均１cmの月刊誌を１年分積み上げたら、高さは平均で何cmになるか？』
答は当然、１２cm。
では次に、
『その雑誌のばらつきは±１mmだった。１年分積み上げたとき、ばらつきは±何mmになるか？』
ここでの「ばらつきが±１mm」の意味ですが、話を単純化して、雑誌の厚さには 9mm と 11mm の２種類しか無くて、
その２種類が毎月ランダムに発行されるものとしましょう。
答の選択肢：１年分のばらつきは・・・
(a) ±１２mm、 ±１mm × 12ヶ月だから。
(b) ±０mm、 ＋と－が打ち消し合って０になるので。
(c) ±１２mm でも ±０mm でもない、中途半端な数。
・・・ちょっと考えてみよう・・・
正解は (c)。
試しに ±1 の乱数を１２個、何回も足し合わせた平均の絶対値をとると 2.7mm くらいの値となります。
（２項係数を用いて詳しく計算すると 2.70703125... といった値になります。）
つまり、ばらつきに関しては単純な足し算ができないということです。
何とか上手い具合に“ばらつきの足し算”ができないものだろうか・・・
そこで編み出されたのが「２乗すれば足し算できる」というアイデアです。
±１mm の乱数が１２個あったとして、それぞれをただ足し合わせるのではなく、２乗して足し合わせてみましょう。
(±1)^2 + (±1)^2 + (±1)^2 + ･･･ ±1)^2 = 12
これなら個々のばらつきがプラスであってもマイナスであっても、合計は確実に+12となるでしょう。
このように、２乗した値によってばらつきを数えたのが、分散という指標だったのです。
・雑誌１冊あたりの分散 ＝ １
・12冊積み重ねた全体の分散 ＝ １２
※12冊のサンプルがたくさんあったときの分散。つまり毎年の雑誌の山をたくさん並べたときのばらつき。
分散は足し算することができる。
この性質を分散の加法性と言うことがあります。
むしろ、足し算ができるような指標を探したところ、うまくいったのが分散であった、というのが実情でしょう。
では、雑誌１年分の分散が12だったとして、実際に雑誌の山の高さは±何mm程度なのでしょうか。
そこで分散の平方根をとって、元の次元に戻したのが標準偏差です。
今の場合、
(標準偏差) = √12 = 3.46mm
といった値となります。
※ 標準偏差 3.46mm は、少し上に書いた平均偏差 2.71mm とは異なります。
※ 標準偏差は偏差２乗和の平均、平均偏差は偏差絶対値の平均です ← しつこいようですが、ここんとこよく間違うので。
■ 分散のエネルギー解釈
・分散とは「バネに溜まったエネルギー」のこと。
・分散とは「足し算できるように」編み出されたばらつきの指標。
この２つが「なぜ分散が２乗なのか」についての主な答になると思います。
さらに穿った見方をするならば、この２つは別々の事柄ではなく、共通の内容を含んでいると解釈することもできるでしょう。
そもそもエネルギーとは何だったかといえば、
「物理的に、時間変化に対して不変な保存量」
のことでした。
何らかのデータ操作について、例えば上に述べた“バネと棒”のような物理的なモデルを考えたなら、
そのモデル上で足し算できる保存量は、自ずとエネルギーに対応付けられるはずです。
こと分散については、「変位の２乗で溜まるエネルギーのような保存量」というイメージで把握できると思います。
この線で標準偏差に解釈を加えるなら、
「各データにエネルギーを等しく分配したとき、１個あたりのデータが変位する大きさ」
ということになります。
■ 最小絶対値法ではダメなのか
以上、分散（と標準偏差）が、いかにばらつきの指標として合理的であるかを見てきました。
ならば、ばらつきの指標は分散だけが唯一無二であり、それ以外にはあり得ないのでしょうか。
実はそうでも無いのです。
分散は、一般的には最も優れた指標なのですが、場合によっては分散以外の指標が有用なこともあります。
先に、分散 ～ 最小２乗法に基づいてデータの当てはまりを考えたのですが、
次に、偏差の値そのもの ～ “最小絶対値法”に基づくデータの当てはまりを考えてみましょう。
１次元の場合、もしデータが２つだけだったなら、
もっとも当てはまりの良い点は、２つのデータの間のどこでも良いことになります。
このように最小絶対値法には、結果が１点に落ち着かないといった不定性があります。
これが最小絶対値法 ～ 平均偏差が流行らない１つの理由です。
データの数をもう少し増やしてみましょう。
もしデータが５つだったなら、最小絶対値法によるもっとも当てはまりの良い点は、中央の３番目の点になります。
平均偏差を基準に考えたなら、
・データが奇数個の場合、順位が中央の点となる。
・データが偶数個の場合、順位が中央に最も近い２点の間のどこでもよい。
つまりこれは「中央値」ということです。
※ 中央値の定義としては、データが偶数個の場合、中央に最も近い２点のちょうど中間（平均）となります。
・分散 ：偏差の２乗の平均値 ： 最小２乗法 ： 最小となる点は平均値
・平均偏差 ：偏差の絶対値の平均値 ： 最小絶対値法 ： 最小となる点は中央値
このように対比させてみれば、決して“分散が優れていて、平均偏差が劣っている”わけでは無いのだと分かるでしょう。
さらに興味深いのは、２次元の場合です。
一般によく用いられている最小２乗法ではなく、最小絶対値法による直線の当てはめも、それはそれでアリなのです。
* 最小絶対値法による回帰分析
>> http://www.orsj.or.jp/~archive/pdf/e_mag/Vol.40_02_261.pdf
* 最小二乗と最小絶対値：ロバスト回帰で外れ値分析
>>https://biolab.sakura.ne.jp/robust-regression.html
* Excelで操る！誤差の絶対値和を最小にする
>> http://godfoot.world.coocan.jp/Abs_Min.htm
最初のリンク先に、ラプラスの考えた最小絶対値に基づく方法が載っていました。
この回帰分析規準（1）を最良のものと考え，そのアルゴリズムを最初に提案したのはLaplace［PiereSimmonMarquisdeLaplace，1749－1827］である。
リンク先の説明を参考に、ラプラスの考えた直線当てはめの方法を試してみました。
・原点を通る直線を考える。
・各データを傾きの順番に並べる Y1/X1 >= Y2/X2 >= Y3/X3 ･･･ Yn/Xn
・X1 + X2 + X3 ･･･ + Xr-1 < Xr + Xr+1 + Xr+2 ･･･ + Xn かつ
X1 + X2 + X3 ･･･ + Xr > Xr+1 + Xr+2 ･･･ + Xn となる r を探す。
つまり、傾きの大きい方から X を合計していって、合計値が残りの合計をちょうど越える番号 r を探す。
・最も当てはまりのよい傾きとして、Yr/Xr を採用する。
これは当てはめ結果の１例です。
青線が最小２乗法、赤線がラプラスの考えた最小絶対値法です。
この例のように極端な外れ値が混じっていた場合、最小２乗法は外れ値に引っ張られるため当てはまりが悪くなります。
中央値が平均値に比べて外れ値に強いように、最小絶対値法には外れ値に対して頑健な性質があります。
回帰分析のあてあまり評価の指標として、RMSE, または MAE という値がよく用いられます。
これらは何なのかというと、分散と平均偏差、あるいは平均値と中央値、という並びに位置づけられるものです。
* 平均平方二乗誤差 = RMSE（Root Mean Square Error）
・RMSE の最小化 = 二乗誤差の最小化 = 誤差に正規分布を仮定した場合の最尤推定
* 平均絶対誤差 = MAE（Mean Absolute Error）
・MAE の最小化 = 絶対誤差の最小化 = 誤差にラプラス分布を仮定した場合の最尤推定
-- 精度評価指標と回帰モデルの評価
>> https://funatsu-lab.github.io/open-course-ware/basic-theory/accuracy-index/
あるいは最近流行りの機械学習に登場する L1正則化(Lasso), L2正則化(Ridge) といったキーワードも、やはりこれらの概念の延長上に位置づけられます。
■ ガウス誤差論より
なぜ２乗なのか、その源流を遡ると、ガウスの誤差論にたどり着きます。
確かに２乗は優れた概念なのだが、ならば２乗を支える根拠・原理は何なのか。２乗以外の概念はあり得ないのか。
私は長らくこの疑問を抱いてきたのですが、元祖ガウスの言葉にようやく答を見たように思います。
今さらですが、（推測）統計学とはデータから帰納的に導かれるものであり、
他の数学のように原理原則から演繹的に導かれるものではありません。
その意味で、２乗が唯一絶対の方法では無かったのです。
ただ、いったん２乗を離れて他のあらゆる方法と比べてみると、改めて２乗こそが最も単純で適していると認めざるを得ません。
これが分散を２乗で数える理由であり、今日に至るまで統計学の根底に流れ続ける基盤だったのです。
- 作者: カール・F.ガウス,Carl Friedrich Gauss,飛田武幸,石川耕春
- 出版社/メーカー: 紀伊國屋書店
- 発売日: 1981/05/01
- メディア: 単行本
- この商品を含むブログを見る
# -*- coding: utf-8 -*- # # Laplaceの絶対値回帰 import numpy as np import matplotlib.pyplot as plt # 直線回帰を行う def fit_line(x, y): mx = np.mean(x) my = np.mean(y) myx = np.mean(y * x) mxx = np.mean(x * x) w0 = (myx - my * mx ) / (mxx - mx ** 2) w1 = my - w0 * mx return ( w0, w1 ) # 適当な直線＋乱数データを用意 N_SAMPLE = 10 ALPHA = 0.0 # 切片は０に限定する BETA = 0.3 EPS = 0.1 X = np.random.rand(N_SAMPLE) # Y = ALPHA + BETA * X + np.random.normal(0, EPS, N_SAMPLE) # 正規乱数 # Y = ALPHA + BETA * X + EPS * (np.random.rand(N_SAMPLE) - 0.5) # 一様乱数 Y = ALPHA + BETA * X + EPS * np.random.standard_cauchy(size=N_SAMPLE) # コーシー分布 # この場合、Laplace回帰はけっこういい感じだ # Y = ALPHA + BETA * X + np.random.laplace(loc=0.0, scale=EPS, size=N_SAMPLE) # ラプラス分布 # Laplace回帰が最尤値になるという。 #### # まずは普通の線形回帰をやってみよう lin, seg = fit_line( X, Y ) print( "線形回帰の傾き,切片 = {:.4f}, {:04f}".format( lin, seg ) ) #### # 以下、Laplaceの考えた方法で回帰をやってみる YperX =Y / X # 個々の点の傾きリスト Slp_index = {} # index->傾きリストを作る for idx in range(N_SAMPLE): Slp_index[idx] = YperX[idx] # 傾きの値で降順にソートする x_list = [] idx_list = [] for key, val in sorted(Slp_index.items(), key=lambda x: -x[1]): # print( "{{ {} : {} }}, X={}".format( key, val, X[key] ) ) idx_list.append(key) x_list.append( X[key] ) # Xの前半の合計値が後半の合計値を越える境界を探す r_index = 0 for r in range( N_SAMPLE+1 ): pre = x_list[ 0 : r ] # 前半 post = x_list[ r : len(x_list)+1 ] # 後半 if sum(pre) >= sum(post): # print( "Pre:", pre, "len=", len(pre) ) # print( "Post:", post, "len=", len(post) ) r_index = r -1 # インデックスなので１引く。 break x_esti = idx_list[ r_index ] # print( "x_ans={}, X={}, YperX={}".format( x_ans, X[x_esti], YperX[x_esti] ) ) lap = YperX[x_esti] print( "Laplace回帰の傾き = {:.4f}".format( lap ) ) # グラフ出力 xlist = np.arange(0, 1, 0.01) # 0.01単位で細かくプロットし、線として出力させる ylist_lsq = [ (lin * x + seg) for x in xlist ] # 線形回帰 ylist_lap = [ (lap * x) for x in xlist ] # Laplaceによる回帰 plt.plot(X, Y, 'o') # 元データ plt.plot(xlist, ylist_lsq, color="blue", label="lsm") # 線形回帰 plt.plot(xlist, ylist_lap, color="red", label="laplace") # Laplaceによる回帰 plt.legend(loc='lower right') plt.show()
サイコロの和と多次元立方体の関係
複数個の一様分布の和は、厳密にはどのような分布に従うのか。
２個のサイコロを同時に振ったとき、最も出やすいのは合計が７のとき。
なぜなら、1+6, 2+5, 3+4, 4+3, 5+2, 6+1 と、６通りの出方があるからです。
それに比べて合計が２となるのは、1+1 の１通り。合計が１２となるのは 6+6 の１通りしかありません。
グラフにすれば、このようなテント型になります。
同じようにして、３個のサイコロの和をグラフにまとめると、こんな形になります。
* 計算サイト・サイコロの和の確率
>> http://www.calc-site.com/probabilities/dice_total
さらにサイコロの数を１０個、２０個と増やしていくと、グラフは滑らかな釣り鐘型曲線へと近づきます。
この釣り鐘型のカーブが「正規分布」なのだ、という話は統計学の教科書に譲るとして、
ここで問題にしたいのは正規分布に至るまでの途中の段階です。
たとえば上のサイコロが３個 ～ 一様分布を３個足し合わせた結果は、厳密にはどのようなカーブなのでしょうか。
あるいは一様分布を７個足し合わせた結果は、正確にはどのような分布なのでしょうか。
先に結果を示します。
以下は、Mathematicaというソフトで一様分布を１～７個まで足し合わせた結果を描いたものです。
（※ただしグラフの高さは正規化されていません。）
一定の美しい規則に従って、釣り鐘型曲線が削り出される様子が見て取れたでしょうか。
Ｎ個の一様分布の和は、Ｎ本の（Ｎ－１）次曲線の組み合わせから成り立ちます。
たとえばＮ＝３（サイコロ３個）の場合、次の３本の２次曲線で、３つの区間をカバーしています。
y = x^2 -- {x = 0～√3/3 の区間}
y = x^2 - 3 ( (x-1(1/√3))^2 ) -- {x = √3/3～2√3/3 の区間}
y = x^2 - 3 ( (x-1(1/√3))^2 ) + 3 ( (x-2(1/√3))^2 ) -- {x = 2√3/3～√3 の区間}
では、どうやってこれらの曲線を求めたのか。
実はこれらの曲線は、Ｎ次元の立方体を、対角線に直行して切った切断面の面積のグラフなのです。
まず２個のサイコロで考えてみましょう。
１個目のサイコロの目を横軸に、２個目のサイコロを縦軸に、全ての出目を図示すると６×６の正方形となります。
この正方形の上で、２つの目の合計が一定となる線は、右下がり45度の直線となります。
たとえば、上の図で青で引いた線は合計が４となる直線。
赤で引いた線は合計が７となる直線です。
そして、これらの合計が一定となる直線の長さは、その合計値が出る場合の数＝確率に比例しています。
そこで、改めて正方形の黄色い対角線を横軸にして見れば（つまり図を45度傾けて見れば）、
正方形がそのままテント型の分布を形作っていることに気付きます。
次に次元を１つ上げて、３個のサイコロの和を考えてみましょう。
３個のサイコロを３次元の縦、横、高さに配置すれば、６×６×６の立方体となります。
この立方体の上で、３つの目の合計が一定となる面は、立方体の対角線に直行する面となります。
立方体の対角線を軸にして、切断面の面積をグラフに描けば、３個のサイコロ（一様分布）の和の分布となります。
立方体の切断面は、図の手前から順に３つのパートに分かれます。
１番手前の青色のパートは、１つの頂点からスタートし、三角錐がどんどん大きくなる過程。
２番目の中間にある黄色のパートは、断面の正三角形から、もう一方の反対向きの正三角形に向けて、６角形の断面が徐々に変形する過程。
３番目の赤色のパートは、１番目と反対に三角錐がどんどん小さくなる過程です。
１番目と３番目のパートが２次曲線となることは、断面積が辺の長さの２乗となることから想像が付くでしょう。
ちょっと悩むのは２番目の、中間にある６角形断面のパートです。
この２枚のスケッチは、６角形断面のパートを取り出して描いたものです。
特に２枚目のスケッチをよく見てください。
断面の６角形とは、拡大する１つの正三角形の３つの頂点から、３つの正三角形を差し引いた残りだと見なせます。
拡大する１つの正三角形と、差し引く３つの正三角形は、実は全く同じ図形を空間的にずらしたものです。
というのも、差し引く３つの正三角形は、元の立方体の隣に位置する、同じ形の立方体の断面だからです。
空間を埋め尽くしている立方体の群れを切断することを想像すれば、
とある１つの立方体の断面積は、同じ形の正三角形を足したり引いたりすることで算出できるのです。
より具体的に、６角形断面の中間パートは次のように定式化できます。
(１番目のパートの正三角形) - ３×(途中から始まった２番目のパートの正三角形）
立方体の対角線の長さは √3 です。
そして、対角線に対する３つのパートの幅はそれぞれ等しく √3/3 ずつです。
なので (途中から始まった…）とは、数字の上では √3/3 だけずらすことを意味します。
以上に基づいて式を完成させたものが、当初に挙げた「一様分布３個」のグラフだったのです。
※３個のサイコロについては、以下がわかりやすいです。
※『ちょっとした「つり鐘型曲線」を描いてみる』
※ https://teenaka.at.webry.info/201704/article_12.html
※『ちょっとした「つり鐘型曲線」と正規分布の関係』
※ https://teenaka.at.webry.info/201704/article_13.html
４個のサイコロ（一様分布）の足し合わせは、４次元立方体を対角線に直行するように切断した断面として表わされます。
４次元立方体を直接描くことはできませんが、断面である３次元図形を想像することは可能です。
４次元立方体の断面は、以下の４つのパートから成り立ちます。
第１パート：正４面体がだんだん大きくなる過程。
第２パート：正４面体の各頂点が削れ、正８面体となる過程。
第３パート：正８面体の４つの面が成長して、正４面体となる過程。
第４パート：正４面体がだんだん小さくなる過程。
これらの過程は、３次元の場合と同じように、最初の正４面体の足し引きによって定式化できます。
第１パート：(正４面体)
第２パート：(正４面体) - ４×(１つずらした正４面体)
第３パート：(正４面体) - ４×(１つずらした正４面体) + ６×(２つずらした正４面体)
第３パート：(正４面体) - ４×(１つずらした正４面体) + ６×(２つずらした正４面体) - ４×(３つずらした正４面体)
式中に登場する正四面体の数、４、６、４ は、二項係数となっています。
二項係数とは、ｎ個の中からｋ個の要素を取り出す組み合わせの数のことです。
なぜ二項係数となるのか。
４次元立方体の各頂点に、４桁の０または１で番号を付けることを考えます。
４つの頂点は４つの次元に対して、それぞれ手前と奥のどちらか（あるいは上下、右左などの）一方に位置付けられます。
最も手前にある頂点 0000 から出発して、1111に向かって辺で結ばれた頂点を線で結んでたどってゆくと、以下のような図形が描かれます。
（このような図を“ハッセ図”と言います。）
この図から、同一切断面に並ぶ頂点の数は「４つの桁の中に何個かの１が入る組み合わせの数」になっていることが分かります。
そしてこの頂点の数は、そのまま足し引きすべき図形（正四面体）の数となっているわけです。
以上で見てきた関係は、そのまま次元数が上がっても変わりません。
・最初のパートはＮ次元の角錐から始まる。グラフの上で、最初の区画は(Ｎ－１)次の曲線となる。
・最初のパートの図形（Ｎ次元の角錐）を、二項係数に従って足し引することによって、２番目以降のパートをつなげる。
（足し引きは、プラスマイナスを交互に繰り返せばよい。）
こうして、Ｎ次元立方体の断面から、Ｎ個の一様分布の和が描き出されます。
※ 最初のパートとなるＮ次元の角錐の体積は、前の記事に記載しました >> d:id:rikunora:20190210
サイコロの数をうんと増やした極限が正規分布になる、ということは、正規分布とは“無限次元立方体の切断面”の形だということになります。
単なるサイコロの足し合わせを、多次元立方体の切断面として捉えるならば、また違った風景が見えてくる。
多次元正三角錐の体積
２次元の正３角形、３次元の正４面体を延長した、４次元空間にある正三角錐の体積はいくつになるか。
さらに延長して、一般にＮ次元の正三角錐の体積はいくつになるか。
４次元空間にある正三角錐のような図形は、正五胞体というのだそうだ >> wikipedia:正五胞体。
ウィキペディアには『超体積: (√5/96)a^4』と書かれているが、これはどうやって計算したのだろうか。
いきなり４次元では想像が付かないので、目に見える２次元、３次元から類推しよう。
まず２次元の正方形、３次元の立方体を考える。
とある１つの頂点から伸びる辺を全て含むように平面で切り取った角錐のことを“コーナーブロック”と呼ぶことにしよう。
“コーナーブロック”とは、上の図で青く塗った部分のことだ。
２次元の場合は直角２等辺三角形、３次元の場合は切り口が正三角形となるような三角錐である。
４次元の場合、コーナーブロックの切り口は正４面体になる。
そして５次元であれば、コーナーブロックの切り口は４次元空間にある正三角錐（正五胞体）となる。
つまり、コーナーブロックの切り口が、知りたかったＮ次元の正三角錐の体積となっているわけだ。
そこで、[コーナーブロックの体積] -> [コーナーブロックの切り口] -> [Ｎ次元の正三角錐] という方針で、
Ｎ次元の正三角錐の体積を計算してみよう。
コーナーブロックが、もとの立方体の何分の１を占めるかを考えてみる。
２次元の正方形の場合には、見ての通り 1/2 だ。
３次元の立方体の場合には、(1/2)×(1/3) = 1/3! となる。
この延長で、４次元のコーナーブロックを考えると (1/2)×(1/3)×(1/4) = 1/4! になる。
なぜ (1/N次元) が掛け合わされるのか。
たとえば３次元空間の角錐は、１つの立方体を３つの方向に、同じ形に切り分けた一部分なので、元の立方体の 1/3 だと分かる。
３次元のコーナーブロックは、２次元のコーナーブロックの底面に対して (1/3) を掛けるので、
結局のところ元の立方体の (1/2) × (1/3) となる。
４次元空間は絵には描けないが、きっと４つの方向に同じ形に切り分けた部分なのだから、
３次元のコーナーブロックの底面体に対して 1/4 を掛けたものになる。
以下、Ｎ次元空間ならコーナーブロックの体積は、元の立方体の 1/N! になる。
次に、コーナーブロックの、立方体の対角線に対する高さを考える。
２次元の（１辺の長さが１の）正方形の対角線の長さは√2、
３次元の（１辺の長さが１の）立方体の対角線の長さは√3、
Ｎ次元の（１辺の長さが１の）超立方体の対角線の長さは√N 。
コーナーブロックの高さは、その対角線の長さの 1/N となる。
２次元であれば、(高さ) = √2/2、
３次元であれば、(高さ) = √3/3 となる。
なぜ、コーナーブロックの高さが、ちょうど対角線の 1/N となっているのか。
その理由は、Ｎ次元立方体が等間隔に空間を埋め尽くすからである。
３次元で考えてみよう。
図の立方体[A]に着目して、１つの頂点から反対側の頂点まで、→a, →b, →c と、３本の辺を伝って移動することを考える。
この →a, →b, →c ３本の辺の長さを、立方体[A] の対角線に投影したとき、３本の長さは等しくなる。
たとえば →b を取り上げてみたとき、立方体[A] の隣にある立方体[B] を想像すれば、
→b とは、ちょうど「立方体[B] にとって、立方体[A] の →a 」と同じ位置にある。
なので、→a と →b の、対角線に投影した長さは等しい。
→c についても同じである。
→c とは、ちょうど「立方体[C] にとって、立方体[A] の →a 」と同じ位置にある。
なので結局のところ、元の立方体[A] の対角線は、→a, →b, →c の投影によって３等分される。
この、「立方体が空間を埋め尽くす」という事情は、次元が幾つ上がっても変わらない（だろう）。
なので、Ｎ次元のコーナーブロックの高さは、いつでも対角線の 1/N となるわけだ。
以上で、コーナーブロックの体積と高さが分かった。
この体積を高さで微分することによって、切り口の面積(体積?!)を知ることができる。
まず馴染み深い(?!)３次元空間でやってみる。
グラフの横軸xは対角線方向の長さ（高さ）、縦軸は切り口の面積S、Sの積分がコーナーブロックの体積Vに相当する。
３次元の場合、切り口の面積は高さの２乗に比例するので、S = a x^2 と置く。（今の時点で a は未知の定数。）
あとは高さxと、体積Vに、上の値を入れて a を計算すればよい。
∫[0～(√3/3)]{ a x^2 }dx = 1/3!
a/3 (√3/3)^3 = 1/6
a = 3√3 / 2
切り口の面積Sは、
a x^2 = (3√3 / 2) (√3/3)^2 = √3/2
これが切り口の正三角形の面積に一致するはずだ。
一辺が１の正三角形の面積は √3/4 なのだが、今考えているコーナブロックの切り口にある正三角形の１辺の長さは√2である。
なので、
S = (√3/4) (√2)^2 = √3/2
確かに一致しているので、どうやらこの方法で良さそうだ。
上と同じことを４次元でやってみよう。
この場合、切り口S には正４面体の体積が現れる。
∫[0～(√4/4)]{ a x^3 }dx = 1/4!
a/4 (1/2)^4 = 1/24
a = 8/3
S = a x^3 = (8/3) (√4/4)^3 = 1/3
１辺が１の正４面体の体積は √2/12。
今求めた切り口S は、一辺が√2の正四面体なので、(√2/12) (√2)^3 = 1/3
確かに一致している。
さらに、５次元でやってみよう。
これで、知りたかった４次元正三角錐の体積が得られる。
∫[0～(√5/5)]{ a x^4 }dx = 1/5!
a/5 (√5/5)^5 = 1/120
a = 25√5/24
S = a x^4 = (25√5/24) (√5/5)^4 = √5/24
これは一辺が√2の４次元正三角錐なので、一辺が１の体積は、
S = (√5/24) / (√2)^4 = √5/96
確かに、Wikipediaに書かれている値と一致した。
以上を一般化して、(N-1)次元の正三角錐の体積が知りたかったなら、
∫[0～(√N/N)]{ a x^(N-1) }dx = 1/N!
a = (√N)^N / (N-1)!
S = a (√N/N)^(N-1) = N^N / {(N-1)! N^(N-1) √N} <-- 一辺が√2 の(N-1)次元の正三角錐
一辺の長さを１にするには、このS を (√2)^(N-1) で割ればよい。
V = S / (√2)^(N-1) = N^N / {(N-1)! (N √2)^(N-1) √N}
正規分布以外の区間推定
ｔ分布を用いた平均値の区間推定には、「正規母集団を仮定する」というただし書きが付いています。
ということは、母集団が正規分布でなかった場合、区間推定は使えないのでしょうか。
母集団の分布は未知であることが多いので、もし正規母集団という条件が絶対だったなら、
区間推定は実際の現場ではほとんど使えない理屈ということになってしまいます。
結論から言えば、たとえ母集団が正規分布で無くても、
母集団から抜き出した標本の平均値が正規分布に従うなら、区間推定の方法を使うことができます。
中心極限定理によれば、ほとんどの分布は（コーシー分布やべき分布といった“平均の無い”分布を除けば）
標本の平均が正規分布に近づくので、区間推定は正規母集団以外にも広く適用可能なのです。
ただし、正規分布から外れるほど大きなサンプルサイズ（データの数）が必要です。
目安は以下の通りです。
・一様分布なら３０個以上
・ベータ分布(a=4,b=2)くらい歪んだ分布なら１００個以上
・指数分布のように一方的に偏った分布なら３００個以上
これらの目安はどうやって得られたのかというと、ひたすらシミュレーションを繰り返した結果です。
以下に掲げる図は、正規分布以外の特定の分布から、何度も繰り返しサンプルを抜き出して区間推定を行った結果です。
図に示す曲線が元になった母集団の分布。
赤でぼんやり描かれているのは、標本を抜き出す操作を１００回繰り返したときの、標本平均のヒストグラムです。
１００回の繰り返しを１０セット行い、ヒストグラムは１０枚重ねて描いてあります。
●正規分布
まずこれが正規分布、サンプルサイズ３０個の場合。
95%信頼区間であれば、100回のうち95回は真の平均値が信頼区間の中に含まれます。逆に言えば、5回は外れます。
データ30個の場合、外れ回数は5.2回。標本平均の標準偏差=0.18。
●一様分布
一様分布、データ30個の場合と、100個の場合。
データ数が増えるほど平均値の推定幅が小さくなる様子が見てとれます。
データ30個の場合、外れ回数は4.5回。標本平均の標準偏差=0.10。
データ100個の場合、外れ回数は5.2回。標本平均の標準偏差=0.06。
データ数30個であっても区間推定は十分有用でした。
●ベータ分布(4,2)
ベータ分布(4,2)、真の平均値は0.666...、歪んだ分布の例として取り上げました。
データ30個の場合、外れ回数は4.2回。標本平均の標準偏差=0.031。
データ100個の場合、外れ回数は5.5回。標本平均の標準偏差=0.018。
データ300個の場合、外れ回数は4.9回。標本平均の標準偏差=0.010。
30個だと結果が不安定（ばらつきが大きい）でしたが、100個あればかなり安定します。
●指数分布
指数分布、真の平均値は1。
データ30個の場合、外れ回数は8.7回。標本平均の標準偏差=0.18。
データ100個の場合、外れ回数は6.6回。標本平均の標準偏差=0.10。
データ300個の場合、外れ回数は5.3回。標本平均の標準偏差=0.06。
元の分布がここまで偏っていたなら、安定した結果には300個のデータが欲しいところです。
●コーシー分布
コーシー分布には“平均値がありません”。なので、区間推定は使えません。
上はデータ30個の場合ですが、確かに見ての通り、区間推定は適用できません。
* コーシー分布とarctan微分 >> [id:rikunora:20171210]
検索したところ、ここに詳しく書かれていました。
* 大人の知らない統計学 ６ 非正規分布の信頼区間
>> https://blogs.yahoo.co.jp/pironotakarabako/63655916.html
以下にシミュレーションのソースを示します。
（分布や数値の付け替えはコメントで行っています。わかりずらい。）
シミュレーションといっても、要は区間推定をひたすら繰り返しているだけです。
# -*- coding: utf-8 -*- """ 区間推定の頑健性、正規母集団以外では、どれほどサンプルがあれば妥当なのか。 """ import numpy as np from scipy import stats import matplotlib.pyplot as plt class IntervalCheck: N_TRIAL_SET = 10 # 試行セット回数 N_TRIAL = 100 # 試行回数 N_SAMPLE = 30 # 100 # 300 # サンプルサイズ＝データ数 up_miss = 0 # 上側に失敗した回数 down_miss = 0 # 下側に失敗した回数 x_min = -2.0 # -1.5 # 0.0 # グラフの範囲 x_max = +2.0 # +1.5 # 1.0 y_scale = 20 # 10 # 分布曲線をｙ方向に拡大して描く # N_TRIAL_SET だけ試行セットを繰り返す def run(self): misses = np.zeros( self.N_TRIAL_SET ) for i in range( self.N_TRIAL_SET ): misses[i] = me.each_run() print( u"信頼区間からの外れ, mean={:.03f}, stdev={:.03f}".format( np.mean(misses), np.std(misses) ) ) plt.xlim( self.x_min, self.x_max ) self.draw_dist() # 分布曲線を描き足す plt.show() return # 分布曲線を描く def draw_dist(self): n = np.linspace( self.x_min, self.x_max, 1000 ) p = [] for i in range(len(n)): x = stats.norm.pdf( x=n[i], loc=0, scale=1 ) # 正規分布 # x = stats.uniform.pdf( n[i], loc=-1, scale=2 ) # 一様分布 # x = stats.beta.pdf( n[i], 4, 2 ) # ベータ分布 # x = stats.expon.pdf( n[i], loc=0, scale=1 ) # 指数分布 # x = stats.cauchy.pdf( n[i] ) # コーシー分布 p.append( self.y_scale * x ) plt.plot( n, p ) return # １回の試行セット def each_run(self): # print( 'sample_size={}'.format(self.N_SAMPLE) ) self.up_miss = self.down_miss = 0 means = np.zeros( self.N_TRIAL ) for i in range(self.N_TRIAL): means[i] = self.trial() miss_ratio = (self.up_miss + self.down_miss) / self.N_TRIAL print( 'miss_ratio={:.02f}, up_miss={}, down_miss={}, '.format( miss_ratio, self.up_miss, self.down_miss ), end="") plt.hist( means, color='r', alpha=0.1, histtype='stepfilled' ) # 半透明で重ね書き print( 'mean_means={:.05f}, std_means={:.05f}'.format( np.mean(means), np.std(means) ) ) return miss_ratio # １回の試行 def trial(self): sample = np.random.normal( 0, 1, self.N_SAMPLE ) # N個の正規分布 (平均, 分散, 出力数) # sample = np.random.rand( self.N_SAMPLE ) * 2 - 1 # N個の一様乱数 # sample = np.random.beta( 4, 2, self.N_SAMPLE ) # β分布、非対称でやってみよう # sample = np.random.exponential( scale=1.0, size=self.N_SAMPLE ) # 指数分布 # sample = np.random.standard_cauchy( size=self.N_SAMPLE ) # コーシー分布は平均が無い real_mean = 0 # 理論上の平均 # real_mean = 4/(4+2) # a/(a+b) ベータ分布の平均はこうなる # real_mean = 1 # 指数分布の平均はscale # ｔ分布による区間推定 mean = np.mean(sample) # 標本平均 t_dist = stats.t( loc = mean, scale = np.sqrt( np.var(sample, ddof=1) / self.N_SAMPLE ), # √(不偏分散/サンプルサイズ) df= self.N_SAMPLE - 1 ) # 自由度 down, up = t_dist.interval(alpha=0.95) # print( '95% interval: {:.3f} < x < {:.3f}'.format(down, up) ) # 区間推定の結果が当たっているかどうか if up < real_mean : self.up_miss += 1 elif down > real_mean : self.down_miss += 1 return mean # 標本平均を返す if __name__ == '__main__': me = IntervalCheck() me.run()
今年はブログを書くヒマがほとんど無かった。来年はもっと時間のゆとりが欲しい。
それでは、よいお年を！
エクセル近似曲線の罠
エクセル近似曲線の「指数近似」「累乗近似」は、いわゆる非線形最小二乗法ではない。
・エクセルで用いているのは、データの対数に直線をあてはめるという方法。
・いわゆる非線形最小二乗法とは、残差（誤差）の２乗和を最小にする方法。
詳しいことは以下のページで尽きているのだが、エクセルの近似曲線は便利だと、
私自身多くの人に勧めている手前、注意を忘れないよう記載しておく。
* Excelグラフ累乗，指数，多項式近似の論文記載の注意(生物科学研究所 井口研究室)
>> https://biolab.sakura.ne.jp/excel-graph.html
* 指数近似
赤い点で描かれているのが元になるデータ、
Y = 0.8 * exp( 0.3 * X ) に、Xに比例する大きさで正規乱数を加えたもの。
下側の青い線がエクセルの指数近似曲線。
上側の緑の線はＲ言語による非線形回帰。
* 累乗近似
赤い点で描かれているのが元になるデータ、
Y = 0.3 * X ^ 1.5 に正規乱数を加えたもの。
直線に近い青い線がエクセルの累乗近似曲線。
よりカーブのきつい、緑の線はＲ言語による非線形回帰。
ぱっと見に、あてはまりが良さそうなのはＲ言語の方に思える。
だからといってエクセルが“間違っている”というわけでもない。
エクセルはＲ言語と異なる方法で線を引いているのである。
Ｒ言語がもともとの回帰式の上で誤差を最小化するのに対し、
エクセルは両辺対数をとった式の上で誤差を最小化する。
* 指数近似
Y = a * X ^ b -- もともとの回帰式
ln(Y) = ln(a) + b ln(X) -- 両辺対数をとった式
* 累乗近似
Y = a * exp(b * X) -- もともとの回帰式
ln(Y) = ln(a) + b * X -- 両辺対数をとった式
もともとの回帰式は曲線だが、両辺対数をとった式は直線となる。
誤差の最小化は、直線の方がずっと簡単だ。
直線に直した上での回帰は以下のようになる。
* 指数近似
* 累乗近似
この直線で求めたパラメータが、エクセル近似曲線のものに一致していることが見て取れる。
※ exp(0.2526)=1.2874, exp(-0.1613)=0.851
エクセル近似曲線は、直線回帰にちょっとオマケを付け足して実現した機能だったのだ。
エクセルの指数近似、累乗近似では、ゼロやマイナスを含むデータに線を引くことができない。
なぜできないかというと、ゼロやマイナスの対数がとれないからである。
また、エクセルには指数回帰を行う"LOGEST"という関数があるが、この関数で求めた数字も、近似曲線で求めた数字と同じである。
（もう１つの、直線回帰のための"LINEST"関数は特に問題無い。）
さらに、エクセルでは近似曲線に対して決定係数R2を計算することができるのだが、
決定係数R2は本来、線形回帰についてのあてはまり指標であり、非線形回帰に用いることはできない。
では、エクセルに出てくる決定係数は何なのかというと、対数をとって直線に直したときの決定係数だったのだ。
なので、この数字をうのみにして、そのまま用いるべきではない。
非線形回帰であてはまりの良さを評価するには、決定係数ではなく、
残差の２乗和、あるいは（回帰分析の）標準誤差を用いるのが妥当だ。
以上の懸念は、指数近似、累乗近似の話であって、他の種類の曲線「対数近似」「多項式近似」では心配しなくてよい。
実際試してみると、対数近似、多項式近似の結果はＲ言語と一致する。
なぜそうなるかというと、
・対数近似，多項式近似 は 線形回帰
・指数近似，累乗近似 は 非線形回帰
だからである。
形の上で、線形＝直線、非線形＝曲線、と思い込んでいる人もいるかもしれないが、
回帰分析の場合、回帰式が推定パラメータについて１次のものが線形、そうでないものが非線形である。
・対数近似の回帰式： y = a ln(X) + b
・多項式近似の回帰式: y = a X^2 + b X + c
いずれの回帰式も（Xではなく）a, b, c については１次式となっている。
で、指数近似，累乗近似がやりたかったとき、結局のところどうすれば良いのか。
・エクセルを使ったときは、エクセルの結果だと明記しておく。
・できればＲ言語、あるいは信頼できる統計用ソフトを使った方が良い。
・実はエクセルの「ソルバー」機能を使えば、非線形回帰も不可能ではない。
でも、そこまでするなら素直にＲ言語使えば良いのではないかな。
* Excelのチャートオプション"近似曲線の追加"機能の評価
>> https://ci.nii.ac.jp/naid/110007025777
以下はＲ言語上での操作です。
### 累乗近似、非線形回帰 # データを用意 x <- c(1:12) y <- c( 1.075513124, 2.104125184, 2.045615008, 3.70051273, 2.221801362, 2.769819828, 6.656385201, 7.037107115, 7.976938733, 8.150842204, 11.70331758, 13.66264136 ) nonline <- nls( y ~ a * x^b, ,start=c(a=1,b=1)) # 非線形回帰 summary(nonline) # 結果の要約 Formula: y ~ a * x^b Parameters: Estimate Std. Error t value Pr(>|t|) a 0.2680 0.1204 2.226 0.0502 . b 1.5587 0.1959 7.957 1.23e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.132 on 10 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 2.96e-06 # プロットしてみる a <- coef(nonline)[1] b <- coef(nonline)[2] plot( y ~ x) curve( a*x^b, add=T, col="red") ### 指数近似、非線形回帰 # データを用意 x <- c(1:12) y <- c( 3.283583359, 2.754235941, 1.674398192, 3.556187361, 0.976339451, 4.652897765, 15.21122786, 12.09124224, 8.109277666, 13.59872153, 23.046961, 24.17609857 ) nonline <- nls( y ~ a * exp(b*x), ,start=c(a=1,b=1)) # 非線形回帰 summary(nonline) # 結果の要約 Formula: y ~ a * exp(b * x) Parameters: Estimate Std. Error t value Pr(>|t|) a 1.53717 0.67319 2.283 0.045521 * b 0.23269 0.04152 5.605 0.000226 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.404 on 10 degrees of freedom Number of iterations to convergence: 15 Achieved convergence tolerance: 3.005e-07 # プロットしてみる a <- coef(nonline)[1] b <- coef(nonline)[2] plot( y ~ x) curve(a * exp(b*x), add=T, col="red")