1 2007/5/14 地球惑星環境学基礎演習 II （田近・小澤・安形）
統計解析（II） 確率密度関数・さまざまな確率分布
（今回の担当：安形＠新領域創成科学研究科自然環境学専攻 [email protected]）0. はじめに
今日は具体的な確率密度関数を学び，それらの関数形および（近似的な）累積関数形を計算する． そしてそれらがどのような用途に用いられるか概観し，あわせて簡単な数表を作成する．1. 確率密度関数と累積確率∼正規分布を例に
ある変量 xiがどのような分布をしているか？ということについては、前回の平均・分散などではまだ 十分な情報でない場合がある。しかし「xiは○○分布である」という一言 1で x iの分布が説明できてしま うことがある。それは、特有の形状をした（そして統計学で古くから研究が進んできた）グラフをもつ 確率分布である。 ここで、xiを連続量とする。そのとき、xiが x と x+dx の間にある確率を f(x)dx とするならば、f(x)は f(x) ≧0 かつ∫
+∞ ∞−
f
(
x
)
dx
= 1
をみたす関数で、これは確率密度関数(Probability Density Function 2 )と呼ばれる。 もっとも有名な確率密度関数は、正規分布(normal distribution)であろう3。これは、 2 2/2 ) (2
1
)
(
σσ
π
m xe
x
f
=
− − （ただし慣例に従って、x の平均を m、分散をσ2で表した）で表現され、N(m, σ2 )と略称される。と くに N(0,1)は標準正規分布という。 1 正確には、どの分布にもパラメタが必要なのでこの「一言」のみでは十分ではないのだが、意味は通じるであろう。 2 筆者は略して PDF と呼んでいたが、Adobe の PDF が登場して以来、確率密度関数を他人に説明するとき話がこんがら がることが多々あった。 3 ガウス分布とも呼ばれる。2 問題(2-1) x, m,σ2が与えられたときに、N(m, σ2 )の x に対応する値を返す関数副プログラム REAL*8 FUNCTION NRMDST (X, M, V)を作成せよ。ここに X,M,V はいずれも REAL*8 型で、それぞれ上記の x, m,σ2の意味である。 さらに、変数 X を-4 から+4 まで 0.01 きざみに増加させ、各 X とそれに対する 標準正規分布 N(0,1)の f(X)をファイル standard-norm.txt に書き出すプログラムを作成せよ。プログ ラムファイル名は normdist1.f とする。また、その結果をグラフに表示せよ。さらに，x=0, 0.31, 0.50, 1. 00, 2.00, 3.144に対する N(0,1)の確率を画面に表示させ，それらがおよそ 0.3989423 0.3802264, 0.3520653, 0.2419707, 0.0539910, 0.0028835 になることを確かめよ．なお，プログラムファイル名は normdist1.f とする。 ＝＝＝＝＝ FUNCTION の作例： c=========================== c Function NRMDST c=========================== real*8 function nrmdst(x, m, v) c=========================== c Declaration of Arguments c implicit none real*8 x, m, v c=========================== c Constants c real*8 PI parameter (PI=3.141592653589793238d0) c=========================== c Function c nrmdst = dexp (-(x - m)**2.0d0 / 2.0d0 / v) $ / $ dsqrt(2.0d0 * PI * v) end ＝＝＝＝＝ここで，これまでに紹介していないテクニックを使っているので解説する：
・ メインプログラム，FUNCTION や SUBROUTINE の冒頭で implicit none を書くと，変数の宣 言が強制される．このため，変数のミスタイプ（これはコンパイルの時にはエラーにならない） によるバグを防げる． ・ 数値のあとに”d0”をつけると，強制的に REAL*8 型実数であるとみなされる56． ・ PARAMETER 文：PARAMETER(名前＝値，名前=値・・・)というフォーマット7．定数を名前で参 照できる8． ・ 実例を出すのは初めてかもしれないが，各行の６文字目が空白でないと，その行は前の行からの 4 これらの x の値自体はそれほど意味を持ってえらんだわけではない． 5
REAL*8 型仮引数をとる FUNCTION や SUBROUTINE に数値定数を与えて呼び出す場合，これが必要になる場合がある． 6
組み込み関数を用いることができるシチュエーションでは，dble を用いて，dble(数値)としてもよい． 7
ここで，「値」には普通は数値そのものが入るが，PARAMETER 文（同じ文でも，既出のものでもよい）の中で定義され た定数を使った式でもよい．例：PARMATER(xdim=9, ydim=8, xydim=xdim*ydim)
8 配列の DIMENSION 宣言で配列サイズとして使えるので便利．複数の配列についてサイズを変えることになったとき も PARAMETER 文のみ変えれば済む．C 言語に親しんでいる人は，機能だけから言えば#define によるマクロ定義に 似ていなくもないと感じるだろう．
3 継続行とみなされる9． ・ 計算式が長くなってしまう場合には，構造をはっきりさせるためにこの行継続とインデント（字 下げ）を併用した表記にすることがある10． ・この継続行のしくみを巧妙に利用した SUBROUTINE 定義表記例がある.例11： c c Subroutine c subroutine radiat $ ( I sun, solar, I dust, co2, O langwv, shrtwv $ )
行継続文字を複数使い分けている．ここで，「I」はその行が SUBROUTINE への INPUT の仮因数で あること，「O」は SUBROUTINE からの出力（結果）を格納する仮引数であることを意味する．これは FORTRAN の仕様では全くないが，コメント的な役割をはたし，あとからプログラムを読むことが容易 になる． ・ 余談：v＝σ2だから dsqrt(2.0d0 * PI * v)は(dsqrt(2.0d0 * PI) * dsqrt(v))と書 いてもいいのだが，そうすると計算負荷の高い dsqrt を二回行なうことになってしまう．しかし それを言うなら，dsqrt(2.0d0*PI)に相当する値も PARAMETER で定数宣言してしまい，たと えばその値を PI2SQR とすると，(PI2SQR * dsqrt(v))とエレガントに書くことも可能である． ただしそのぶんプログラムは直感的ではなくなる． 9 かつての FORTRAN には一行の長さに厳しい制限があった．今でも，72 文字以上の長さの行は好ましくない．しかし， この文継続を用いれば，継続した行が合計 72 字を越えても問題ない． 10 ただし個人の好みがあるのは確かである．長い継続行をつくらず，細かい途中計算をすべて別々の変数に代入し，最 後にまとめて計算するという方法のほうが好まれる場合がある．とくに式が極度に複雑な場合には後者のほうがデバグ しやすい． 11 実在のプログラムから抜粋したものだが，仮引数名は変えてある．
4 ＝＝＝＝＝ 問題(2-2)：x が与えられたとき N(0,1)の f(x)を与える関数副プログラム REAL*8 FUNCTION STDNRM(x) ( x は REAL*8 型)を作成せよ．なお，FUNCTION 内部から上記の NRMDST を呼び出してもよい． 作例：主要部分のみ書くと， implicit none real*8 x real*8 nrmdst stdnrm = nrmdst(x, 0.0d0, 1.0d0) ========================= さて，確率密度関数は，「ある確率密度関数に従う xiの値が x となる確率」をあらわすわけだが，実用 的には（これは本日および次回の内容でもある）， ・ ある確率密度関数に従う x の値が a 以上(ないし a 以下，あるいは|xi|≧a)になる確率は何か や ・ x がある確率密度関数に従うとき，x の値が a 以上(ないし a 以下，あるいは|x|≧a)になる確率が○ ○%になるような a の値は何か といった問題によく使われる．このためには，次の分布関数 F(t)が使われる12：
∫
−∞=
ｔ(
)
)
(
t
f
x
dx
F
これは確率密度関数 f(x)に従う確率変数が x 以下となる確率をあらわす13．明らかに以下の式が成り立 つ： F(x)≧0，(
)
≥
0
dx
x
dF
，lim
(
)
=
0
−∞ →F
x
x ，lim
(
)
=
1
∞ →F
x
x また，正規分布のような平均値 m の周りに左右対称な確率密度関数ならば， F(m+x)+F(m-x)=1, F(m)=0.5 例題(2-1) [a]N(0,1)について，台形公式（刻み幅 0.01）による f(x)の数値積分を用いて，0.00≦x≦4.00 （ただし x は 0.01 刻みで動かすものとする）について F(x)を求めて画面とファイル‘erf1.txt’に出 力する（ただし）プログラム PROGRAM ERF1(ファイル名は erf1.f)を作成せよ．上記の関数副プロ グラム STDNRM を用いてよい． 12 混乱を避けるために F(x)ではなくて F(t)と表記したが実際には F(x)と書かれることは多い． 13 グラフ y=f(x)上で考えるならば,ｘ軸，直線 y=x で囲まれる部分の面積である．5 作例： program ERF1 c c Variables c implicit none real*8 stdnrm
real*8 x, minx, maxx, step real*8 nd, sum, prev integer i, ifirst, ilast c
c Constants c
parameter (maxx=4.0d0, minx=0.01d0, step=0.01d0) c c Main c c c Open File c open(1, FILE='erf1.txt') c c Initialize c
ifirst = int(minx / step) ilast = int(maxx / step) sum = 0.5
prev = stdnrm(0.0d0) c
c write first value (for x=0.0) c write(6,1000)0.0, sum write(1,1000)0.0, sum c c Loop c do i = ifirst, ilast x = dble(i) * step nd = stdnrm(x) sum = sum + $ ( step * ( nd + prev ) / 2 ) write(6, 1000)x, sum write(1, 1000)x, sum c c nd (f(x) ) --> prev c prev = nd end do 1000 format(F6.2, F12.8) c c File Close c close(1) stop end
6
ここでもいくつかプログラム上のテクニックを使っている：
・ x が変化する最初と最後の値(maxx, minx)およびその動く幅(step)を PARAMETER 文にて定数とし て宣言している．
・ ループ変数 i の最初と最後の数値をこの maxx, minx, step から計算している．
・ 台形積分の公式をもちいる際，各ループにおいて f(x)と f(x-Δx)を両方計算しなおすのではなく， 後者は前回のループにおける f(x)に等しいことを利用してそれを prev という変数に入れて参照し ている． ・ これは復習に近いが，一つの FORMAT 文を複数の WRITE 文で共有することができる． 問題(2-3) 上記により作成されたファイル erf1.txt を http://hydro.iis.u-tokyo.ac.jp/~agata/doc/2007KisoEnshu/20070514.html に置いた．この内容から，N(0,1)における F(1.0), F(2.0), F(3.0)の概数を求めよ．また， F(x)=0.995， F(x)=0.975 となる x の値を推定せよ14．これらの値は何を意味するか． =========
問題(2-4a) 上記のプログラム erf1.f をコピーして erf2.f を作成せよ．そして，主プログラムを 改造して x<0 の場合にも対応するようにせよ．結果は-4≦x≦4 の範囲で x を 0.01 刻みで動かして，f(x), F(x)の値をそれぞれ配列 NDST, ERFA(添え字の意味は，ERFA(1)=F(0.01), ERFA(2)=F(0.02)…とする) に格納したのち，x と F(x)の組を例題 2-1 と同様に画面及びファイル nd-and-erf.dat に出力せよ． 問題(2-4b) ここまでに得られた，N(0,1)の f(x)と F(x)を同一平面上にグラフ化せよ このようになっただろうか？ 参考コードが次ページにあるが，自力で取り組んでみてほしい． 14 台形公式というなんとも誤差の大きそうな近似を用いているが，F(x)の例では相対誤差は 10-4程度で済んでいる．
7 参考コード（erf2.f 主プログラムの抜粋）
c
c Initialize c
ifirst = int(minx / step) ilast = int(maxx / step) sum = 0.5 prev = stdnrm(0.0d0) c c Loop (for x>0) c do i = ifirst, ilast x = dble(i) * step nd = stdnrm(x) sum = sum + $ ( step * ( nd + prev ) / 2 ) c
c Store (tentative) results into arrays c ndist(i) = nd erfa(i) = sum c c nd (f(x) ) --> prev c prev = nd end do c Calculation Finished c c Writing... c c c Open File c open(1, FILE='nd-and-erf.dat') c (x:Negative value) do i=ilast, ifirst, -1 x = - step * dble(i)
write(6,1000)x, ndist(i), 1.0d0 - erfa(i) write(1,1000)x, ndist(i), 1.0d0 - erfa(i) end do c (x:zero) x = 0.0d0 nd = stdnrm(x) write(6,1000)x, nd, 0.5 write(1,1000)x, nd, 0.5 c (x:Positive value) do i=ifirst, ilast x = step * dble(i)
write(6,1000)x, ndist(i), erfa(i) write(1,1000)x, ndist(i), erfa(i) end do 1000 format(F6.2, 2F12.8) c c File Close c close(1)
8 近似計算：F(x)を任意の x について求めようと思ったらどうすればよいだろうか？実際には解析的な方 法も，誤差がまったくない数値解法もないので，近似計算にて行なうしかない．ここでは Hastings の近 似解法を紹介する： Hasting の４次近似式：x≧0 のとき、N(0,1)の F(x)＝1-0.5/w ただし，w=(a0 + a1 x + a2x 2 + a3x 3 + a4x 4 )4で ある 。ここに、a0=1.0 , a1= 0.196854 , a2= 0 .115194, a3= 0 .000344 , a4= 0.019527 とする。 x<0 の ときは F(x)=1-F(| x |)から求めればよい。 例題 2-2 Hastings の４次近似式を用いて N(0,1)の F(x)を任意の x について求める関数 FUNCTION ERFH4(x)を作成せよ。 解１： c========================================================= c c Function ERFH4 c N(0,1)の F(x)の計算（Hastings の 4 次近似）． c
real*8 function erfh4(x) c---
c Declaration of Arguments and Functions c
implicit none real*8 x, x1 real*8 w
real*8 a0, a1, a2, a3, a4
parameter (a0 = 1.0 , a1 = 0.196854 , $ a2 = 0.115194, a3 = 0.000344 , a4 = 0.019527) c--- c Function c x1 = dabs(x) w = a0 + a1 * x1 + a2 * x1 ** 2 + a3 * x1 ** 3 + a4 * x1 ** 4 w = w ** 4.0 if x >= 0 then erfh4 = 1.0 - 0.5 / w else erfh4 = 0.5 / w end if end これでも悪くはないが、実は二つまだ改良できる。ひとつは多項式の計算法。もう一つはデータの大小 によって符号が異なる値の扱い方。組み込み数学関数の sign を使えば if 文を全く使わずに式を書け る。 c========================================================= c c Function ERFH4 c N(0,1)の F(x)の計算（Hastings の 4 次近似）． 改良版
9
c
real*8 function erfh4(x) c---
c Declaration of Arguments and Functions c
implicit none real*8 x real*8 w
real*8 a0, a1, a2, a3, a4
parameter (a0 = 1.0 , a1 = 0.196854 , $ a2 = 0.115194, a3 = 0.000344 , a4 = 0.019527) c--- c Function c c c 多項式はこのようにして計算する c
w = a0 + x * (a1 + x * (a2 + x * (a3 + x * a4))) w = w ** 4.0 c c sign 関数を使った x の正負での計算振り分け c erfh4 = 0.5 + sign(0.5 - 0.5 / w, x) end 問 題 (2-5) こ れ ま で に 作 成 し た 問 題 2.4 の プ ロ グ ラ ム お よ び FUNCTION erfh4 を 用 い て 、 x=-4.00,-3.99,…,3.99,4.00 までの x, f(x), （問題 2.4 の）F(x)、（Hastings の４次近似の）F(x)をファイル nd-and-erf-Hastings4.dat に書き出せ。その結果をグラフにプロットし、二種類の F(x)のグラフ がどれほど重なるか確かめよ。
10
2. 相関から回帰へ
回帰直線：前回の宿題で扱った三つの X,Y データセットのプロット（散布図：Scatter Diagram）は以下 のようになる： (A vs B : r≒0.98143) (A vs C : r ≒0.71645） (A vs D : r ≒0.06311) AvsB や AvsC の散布図を見ると、散らばった点が一本の 直線で代表される、言い換えれば、x とｙの関係をもっと もよく表す直線が定義できそうな気がしてくる。これは 次のように予測偏差を最小にするという定義づけのもと では決定可能である。一般にこれを回帰直線という。 変数ｘiから yiを予測するための回帰直線を y=ax+b とし よう。このとき、xiに対する yiの予測値を yi’（予測のズ レは yi’− yiとなる）とすると、次の二つの条件をともに 満たす直線が「最良」だと考えられるだろう： Σ(yi’−yi)=0 (偏差の平均が０：全体として偏りがない) Σ(yi’−yi) 2が最小 (偏差の平均が０：全体として偏りがない) そこで、Σ(yi’−yi) 2＝Σ(ax i+byi’−yi) 2 が最小となるように a と b を決めればよい。展開し、a および b で偏微分したのちに 0 とおいて、結 局最後はこの式を得る： a（傾き）については 2 2 2
(
)
)
y
)(
(
x
x
y
x
x
x
n
x
y
x
n
y
x
a
i i i i i i−
−
−
=
⋅
−
⋅
⋅
−
=
∑
∑
∑
∑
となり、b（切片）については次の関係式から求められる：y
−
y
=
a
(
x
−
x
)
（これは回帰直線が x,y の平均値となる点を通ることを意味する）x
iy
i’y
i’y
i’―y
i11
問題 2-6：INTEGER 型の N,REAL*8 型の配列 X(1000),Y(1000)を受け取ったとき、相関係数 R, 回 帰直線の傾き A、同切片 B を求める SUBROUTINE REGRES（N,X,Y,R,A,B）を作成せよ。それを用い て、前回宿題の corr.txt にある A と B の間の回帰直線を求めよ。同様に、A と B、A と D について もおこなえ。プログラム名は regres1.f とせよ。
3. 相関・回帰の怖さ
技術的な怖さ： 上に挙げたのは１変数のみで他の変数（被説明変数）を予測しようとする回帰であり、それは１次式 に基づいたもの（直線回帰）であった。このような回帰は次のような怖さを持っている。 ・ むしろ曲線のほうが回帰しやすい場合には、直線の回帰はうまくいかないか、変量間の関係を見 逃しがちになる。 ¾ ただしお互いの変量の log をとったり、変数変換を行なうことにより、直線回帰に帰着でき ることもある。 ・ 本質的に単調増加ないし単調減少でない関係のとき（たとえば、クズネッツ曲線のように馬蹄形 のグラフになる場合）、一つの y に対して対応する x が複数あることになり、回帰の考え方にそぐ わない ・ 相関関係は因果関係ではない。 ¾ これにはいろいろなパターン（落とし穴）がある。 ¾ １．全くの偶然で相関があるようなデータが出てくることがある。 異なる要因がドライブされる二つの変量があり、その要因がたまたまデータ取得期間に のみ同一の傾向で変動した など ¾ ２．同じ要因で、しかし独立に変動する二つの変量 ・ 少数の点が飛びぬけた x,y の組の値をとるとき、見かけ上相関係数が高くなる。相関係数は変量間 の関係についてアタリをつけるには向いているが、散布図も書いてみて慎重に検討する必要があ る。12 ・ 有名な例15： グラフ中の数字は年度。年度ご とに 二つの変量 （統計資料や proxy data などから取得）を求め、 その関係を求めた。 問題 2-7：二つの変量の年度変 化をグラフから読み取り，それ をあ らわすテキ ストファイル carib.dat を作成することに より，corr.f を用いて両者の 相関係数を求めよ． ＝＝＝＝＝ レポート課題：次回まで．＜１＞＜２＞は A4５枚以内．＜３＞＜４＞は無制限．また，電子ファイル を [email protected] まで送付．文書フォーマットは自由だが，PDF ないし DOC が望ましい． ＜１＞ 地球惑星環境学の課題の中から，相関があると思われるデータの組を３組選び，それらの記 述統計（平均，不偏分散，標準偏差，必要があれば中央値，最頻値，歪度，尖度）を記載せ よ．そして両者の相関係数および回帰直線を求めよ．なお，データの集め方は自由だが，か ならずデータソースを示すこと．また，計算に使ったプログラムおよびデータファイルと， 結果のログを添付すること． ＜２＞ ＜１＞の結果は地球惑星環境学の観点からみて何を意味しているか．３組それぞれ，根拠と ともに述べよ． ＜３＞ やや上級者編１：今週中にメイルで指示します． ＜４＞ やや上級者編２：今週中にメイルで指示します． （安形注：３，４は第三回レジュメの末尾参照） 15 http://www.venganza.org/about/open-letter/より 何のグラフだろう？
