今日のテーマ: 相関
今回の目標は、松山市の過去134年間における各月の気温の相関係数を求め、またそれらの間の関係を「あてはめ」た回帰直線の方程式を求めること。
最初の目標として、「過去134年間の5月と12月の月平均気温の相関係数を求める」ことを考えよう。 今日の演習でも、第4回 で作成した、平均・分散・標準偏差を計算するプログラム (enshu03c.f90 という名前だったか) の 適切な場所に、以下に示す「部品」をうまく追加する ことで目的のプログラムが完成するようになっている。 ただし、どこに追加するのが適切かは 各自で考える こと。
また、講義でも指摘したように、相関・相関係数を求めるための手順は、線形最小二乗法で直線に「あてはめ」る手順と非常によく似ている。 その反面、変数の意味の違いによる混乱が起こりやすいので、「最小二乗法」のときに使ったプログラムからのコピペはしないのが無難 だと思う。
x と y の積の平均を求める
必要な手順 (の要点だけ) を以下に示す。
「5月の気温と12月の気温の相関係数」を求めるのに必要な「器」を用意する。 例えば「気温の積の総和」用に seki_souwa、「気温の積の平均」用に seki_heikin という名前の 実数 を宣言する。
r e a l : : s e k i _ s o u w a , s e k i _ h e i k i n 「5月の気温と12月の気温の積の総和」を入れる「器」の中身をあらかじめ「空っぽ」にしておく。
s e k i _ s o u w a = 0 . 0 「気温の積の総和」を入れる「器」にどんどん値を足していく。 各年の各月の気温を読むたびに、seki_souwa という器に足し込んでいる。
d o r e a d ( * , ' ( i 4 , 1 2 f 5 . 1 ) ' , i o s t a t = j o u t a i ) s e i r e k i , k i o n i f ( j o u t a i / = 0 ) e x i t n e n s u u = n e n s u u + 1 s e k i _ s o u w a = s e k i _ s o u w a + k i o n ( 5 ) * k i o n ( 1 2 ) e n d d o 「気温の積の平均」を計算する。 ここでも当然、/nensuu ではなく /real(nensuu) とする 必要がある。
s e k i _ h e i k i n = s e k i _ s o u w a / r e a l ( n e n s u u )
相関係数を求める
ここまでくれば後は簡単、のはず。 必要な手順 (の要点だけ) を以下に示す。
必要な「器」を用意する。 例えば soukan という名前の 実数 を宣言する。
r e a l : : s o u k a n 相関係数 soukan の値を計算し、表示させる。 講義のノートを参照して、適切な計算式を書き入れよ。 この際、5月と12月の標準偏差で2回割り算する手順を /kion_tsuki_hensa(5)*kion_tsuki_hensa(12) と書いてはいけない ことにも注意すること。 また相関係数の値を出力するにあたり、小数点以下を3桁まで表示するようにしてある。
s o u k a n = w r i t e ( * , ' ( f 7 . 3 ) ' ) s o u k a n
これらを参考にして、過去134年間の5月と12月の月平均気温の相関係数を求める手順がどう書けるかを考えよ。 そして 第4回 で作成した、平均・分散・標準偏差を計算するプログラム (enshu03c.f90 だったか) にその機能を追加したプログラムを作成せよ [みほん]。
gfortran に怒られないプログラムができたら、例の如く
$ cat matsuyama-kion.txt | ./a.outと実行して計算結果を確認せよ。 正しくプログラミングできていたら、プログラムの出力 (の最後の1行) は以下のようになっているはずである。
0.479
また、実際に両者の相関の程度をグラフにして眺めてみよう。 端末から Gnuplot を起動し、以下のように入力してみよう (gnuplot> は Gnuplot のプロンプト)。
gnuplot> plot "matsuyama-kion.txt" using 6:13とすると、5月の平均気温を横軸に、12月の平均気温を縦軸にとったグラフが表示されているはずだ。 相関係数の値 (符号、および絶対値の大きさ) と、実際にグラフに表示されているデータの相関の程度を確認してみよ。
余力のある人向けのおまけ: 回帰直線を求める
今度は、上で表示したグラフの回帰直線の方程式や、「あてはめ」の決定係数を求める機能を追加してみよう。 必要な手順 (の要点だけ) を以下に示す。 ただし、すでに相関係数や各月の気温の平均・分散・標準偏差が全て求まっている ので、新たに計算しなければいけないものは何もない。
必要な「器」を用意する。 例えば、直線の傾きを kaiki_a、y切片を kaiki_b、および「あてはめ」の決定係数を r2 に入れることにする。
r e a l : : k a i k i _ a , k a i k i _ b , r 2 「あてはめ」の決定係数、および回帰直線の傾きと y 切片を計算し、表示させる。 講義のノートを参照して、適切な計算式を書き入れよ。 なお当然ながら、決定係数や回帰直線の傾きを求めるには、既に求まっている相関係数の値をできる限り利用した式で計算 するのが簡単。 さらにその際、何が x で、何が y かに十分注意 すること。 ここでは 5月の平均気温を横軸に、12月の平均気温を縦軸に とったグラフを描く場合を例にとっているので、そのつもりで x と y に相当するものを選ぶこと (さもなければ、「式も答も合っているのに、グラフを描かせてみたら全然違う」ことになる)。 また計算結果を出力するにあたり、相関係数と決定係数は小数点以下を3桁まで表示するようにしてある。
w r i t e ( * , ' ( i 4 , 1 2 f 6 . 2 ) ' ) n e n s u u , k i o n _ t s u k i _ h e i k i n w r i t e ( * , ' ( i 4 , 1 2 f 6 . 2 ) ' ) n e n s u u , k i o n _ t s u k i _ b u n s a n w r i t e ( * , ' ( i 4 , 1 2 f 6 . 2 ) ' ) n e n s u u , k i o n _ t s u k i _ h e n s a s e k i _ h e i k i n = s e k i _ s o u w a / r e a l ( n e n s u u ) s o u k a n = r 2 = k a i k i _ a = k a i k i _ b = w r i t e ( * , ' ( 2 f 7 . 3 , f 6 . 3 , f 6 . 1 ) ' ) s o u k a n , r 2 , k a i k i _ a , k a i k i _ b
これらを参考にして、過去134年間の5月と12月の月平均気温の相関係数と回帰直線の方程式を求める手順がどう書けるかを考えよ。 そして、今回の演習で先に作ったプログラムにその機能を追加したプログラムを作成せよ [みほん]。
gfortran に怒られないプログラムができたら、例の如く
$ cat matsuyama-kion.txt | ./a.outと実行して計算結果を確認せよ。 正しくプログラミングできていたら、プログラムの出力 (の最後の4行) は以下のようになっているはずである。
134 5.26 5.59 8.64 13.72 17.97 21.91 26.24 27.05 23.48 17.69 12.51 7.68 <- 平均 134 1.51 2.10 2.00 1.52 1.33 0.98 1.37 1.09 1.74 1.53 1.60 1.53 <- 分散 134 1.23 1.45 1.42 1.23 1.15 0.99 1.17 1.04 1.32 1.24 1.27 1.24 <- 標準偏差 0.479 0.230 0.513 -1.5 <- 相関係数、決定係数、回帰直線の傾きと y 切片
さらに、直線への「あてはまり」具合をグラフで確認してみよう。 再び Gnuplot を起動し、以下のように入力してみよ。
gnuplot> plot "matsuyama-kion.txt" using 6:13, 0.513*x-1.5このグラフで、回帰直線と実際のデータの位置関係がもっともらしいものになっているか どうか、また 回帰直線への「あてはめ」の程度が、相関係数 (あるいは決定係数) の大きさから予想されるものと合っているか どうかを改めて確認せよ。
相関係数を求める (その2)
さらに今度は、5月の平均気温と、全ての月の平均気温との相関係数や回帰直線の傾き・y切片・決定係数をいっぺんに求めるプログラムを作ってみよう [みほん]。
必要な手順 (の要点の一部) を以下に示す。
必要な「器」を定義する。 「気温の積の平均」を求めるのに必要な変数、および相関係数などなどを、12個の実数を入れる配列 として宣言する。
r e a l , d i m e n s i o n ( 1 2 ) : : s e k i _ s o u w a , s e k i _ h e i k i n , s o u k a n r e a l , d i m e n s i o n ( 1 2 ) : : k a i k i _ a , k a i k i _ b , r 2 各月の「平均気温と5月の平均気温の積の平均」を計算する。 やるべきことは要するに
からs e k i _ s o u w a ( 1 ) = 0 . 0
だったり、またs e k i _ s o u w a ( 1 2 ) = 0 . 0
からs e k i _ s o u w a ( 1 ) = s e k i _ s o u w a ( 1 ) + k i o n ( 5 ) * k i o n ( 1 )
だったり、さらにはs e k i _ s o u w a ( 1 2 ) = s e k i _ s o u w a ( 1 2 ) + k i o n ( 5 ) * k i o n ( 1 2 )
からs e k i _ h e i k i n ( 1 ) = s e k i _ s o u w a ( 1 ) / r e a l ( n e n s u u )
という訳なのだが、まさかこの期に及んで12個の式を長々と書き並べるなんてことのないように。s e k i _ h e i k i n ( 1 2 ) = s e k i _ s o u w a ( 1 2 ) / r e a l ( n e n s u u ) 各月の「平均気温と5月の平均気温との相関係数」および「平均気温と5月の平均気温との回帰直線の傾きとy切片」のを計算する。 5月と各月の相関係数などを計算するには、気温の平均や標準偏差が12ヶ月分全て必要 であるから、この部分は、各月と5月の気温の積の平均だけでなく、各月の気温の平均や標準偏差を全て求める部分より後に追加 するのがよい。
d o i = 1 , 1 2 s o u k a n ( i ) = r 2 ( i ) = k a i k i _ a ( i ) = k a i k i _ b ( i ) = e n d d o 最後に結果を表示する。 この部分は以下のようにするのがよい。
w r i t e ( * , ' ( i 4 , 1 2 f 6 . 2 ) ' ) n e n s u u , k i o n _ t s u k i _ h e i k i n w r i t e ( * , ' ( i 4 , 1 2 f 6 . 2 ) ' ) n e n s u u , k i o n _ t s u k i _ b u n s a n w r i t e ( * , ' ( i 4 , 1 2 f 6 . 2 ) ' ) n e n s u u , k i o n _ t s u k i _ h e n s a w r i t e ( * , ' ( 4 x , 1 2 f 6 . 3 ) ' ) s o u k a n w r i t e ( * , ' ( 4 x , 1 2 f 6 . 3 ) ' ) r 2 w r i t e ( * , ' ( 4 x , 1 2 f 6 . 3 ) ' ) k a i k i _ a w r i t e ( * , ' ( 4 x , 1 2 f 6 . 1 ) ' ) k a i k i _ b
gfortran に怒られないプログラムができたら、例の如く
$ cat matsuyama-kion.txt | ./a.outと実行して計算結果を確認せよ。 正しくプログラミングできていたら、プログラムの出力 (の最後の7行) は以下のようになっているはずである。
134 5.26 5.59 8.64 13.72 17.97 21.91 26.24 27.05 23.48 17.69 12.51 7.68 <- 平均
134 1.51 2.10 2.00 1.52 1.33 0.98 1.37 1.09 1.74 1.53 1.60 1.53 <- 分散
134 1.23 1.45 1.42 1.23 1.15 0.99 1.17 1.04 1.32 1.24 1.27 1.24 <- 標準偏差
0.410 0.585 0.651 0.646 1.000 0.640 0.475 0.662 0.561 0.710 0.636 0.479 <- 相関係数
0.168 0.343 0.424 0.417 1.000 0.410 0.225 0.438 0.315 0.504 0.405 0.230 <- 決定係数
0.437 0.735 0.799 0.690 1.000 0.549 0.482 0.599 0.642 0.761 0.698 0.513 <- 回帰直線の傾き a
-2.6 -7.6 -5.7 1.3 0.0 12.1 17.6 16.3 11.9 4.0 -0.0 -1.5 <- 回帰直線の y 切片 b
正しく動いていることが確認できたら、本日の演習の出席確認のため、このプログラムを 末尾のフォーム 経由で提出せよ。
また、Gnuplot を使って先と同様に、横軸に5月の平均気温、縦軸にある月の平均気温をとったグラフを描いてみよ。 データの相関の度合が、相関係数の値から予想されるものと合っているかどうかを確認せよ。
余力のある人向けのおまけ: 2次元配列・2重ループ
ここまでの課題では、5月の気温と各12ヶ月の気温との相関をとってきた。 その際、5月の気温とi月の気温との相関 を入れておくための変数を、12個の実数を入れる配列 soukan として宣言し、i月の気温との相関係数を 配列 soukan の i 番目の「番地」に入れる ようにしていたのであった。 では、もし全ての月の間の相関、具体的には i月の気温とj月の気温の相関係数を全部 求めてみたくなった場合には、どのような方法を用いればよいだろうか。
Fortran 90 に限らずコンピューター言語のほとんどは複数の「添え字」を使った配列を用いることができる。 具体的には、2つの「添え字」を使った配列を2次元配列 といい、kion(12) のように1つの「添え字」を使った配列を1次元配列という。 以下の例は、変数 seki_souwa、seki_heikin、soukan を2次元配列として宣言するためのものである。 この場合、これら3つの配列は、1つめの「番地」は1から12まで、2つめの「番地」も1から12までの値をとり、合計で 12×12=144 個の実数 を入れることができる。
| r | e | a | l | , | d | i | m | e | n | s | i | o | n | ( | 1 | 2 | , | 1 | 2 | ) | : | : | s | e | k | i | _ | s | o | u | w | a | , | s | e | k | i | _ | h | e | i | k | i | n | , | s | o | u | k | a | n |
余力のある人は、このプログラム (コピペ可) を参考に、全ての月の間で平均気温の相関係数を求めてみるプログラムを作成してみよ。 ある月とそれ以外の月の組み合わせをとった \({}_{12}\mathrm{C}_{2}=\frac{12\times11}{2}=66\) 通りのうち、何月と何月の平均気温の間が最も相関がよいか。
ちなみに、このプログラムの出力 (の最後の13行) は以下のようになるはず。 この「表」形式に書かれている \(12\times12=144\) 個の相関係数の結果からも、
- i月の気温とi月の気温の相関係数は1
- i月の気温とj月の気温の相関係数は、j月の気温とi月の気温の相関係数と同じ
1 2 3 4 5 6 7 8 9 10 11 12
1 1.000 0.468 0.402 0.370 0.410 0.244 0.120 0.280 0.292 0.361 0.371 0.221
2 0.468 1.000 0.580 0.430 0.585 0.425 0.256 0.463 0.386 0.462 0.380 0.357
3 0.402 0.580 1.000 0.518 0.651 0.518 0.392 0.535 0.465 0.553 0.444 0.412
4 0.370 0.430 0.518 1.000 0.646 0.488 0.428 0.513 0.467 0.568 0.371 0.333
5 0.410 0.585 0.651 0.646 1.000 0.640 0.475 0.662 0.561 0.710 0.636 0.479
6 0.244 0.425 0.518 0.488 0.640 1.000 0.496 0.480 0.472 0.561 0.530 0.369
7 0.120 0.256 0.392 0.428 0.475 0.496 1.000 0.575 0.456 0.459 0.287 0.233
8 0.280 0.463 0.535 0.513 0.662 0.480 0.575 1.000 0.616 0.611 0.396 0.235
9 0.292 0.386 0.465 0.467 0.561 0.472 0.456 0.616 1.000 0.656 0.473 0.249
10 0.361 0.462 0.553 0.568 0.710 0.561 0.459 0.611 0.656 1.000 0.551 0.397
11 0.371 0.380 0.444 0.371 0.636 0.530 0.287 0.396 0.473 0.551 1.000 0.492
12 0.221 0.357 0.412 0.333 0.479 0.369 0.233 0.235 0.249 0.397 0.492 1.000
さらに、上で挙げたプログラムでは、2次元配列を操作している部分は、do〜end do のループの中にもう1つの do〜end do のループを含んでいる、という2重構造になっていることにも注意してほしい。 この例のように、複数のループを「入れ子」にして用いる ことも可能であり、2つのループが「入れ子」になったものは2重ループと呼ばれる。 また、このプログラムの出力を注意して眺めることにより、2重ループの中の2つのカウンタ変数 (i と j) がどのような順番で変化しているかを確認せよ。
出席確認用プログラム提出フォーム
Moodle のコースページ 以下の 課題ファイル提出 により、各自で作成した (enshu05c.f90 に相当する) 過去134年間における各月の平均気温と5月の平均気温の相関係数を求める Fortran 90 プログラムの最終版を提出せよ。