計測工学 I 第 4 回 Excel による相関係数の計算
今日の内容 第 4 回 Excel による相関係数の計算 シラバスより以下の部分を扱います 第 4 回 Excel による相関係数の計算 相関の考え方を導入する 相関が高いということ 負の相関があるという概念を学び 時系列データの位相差 相互相関関数 相関係数 自己相関関数について学ぶ 二つのデータが伴って変わる時に どんな関係が成り立つかを求めるのが相関係数です
相関とは? 二つのものが密接にかかわり合っていること - Mac 国語辞典 因果関係を調べる 原因がある 結果につながる : 因果律 自然科学は因果律の塊とも言えるかも知れない A という現象が B という現象に影響する 一定の遅れを伴って 両者が同様の傾向を示す 相関するものの例 身長と体重 標高と気温 試験前の勉強時間と得点 太陽の高さと気温 など
相関を調べる 自然科学では 未知の現象の 相関 を調べるものが多い HbA1C( ヘモグロビン A1C) は過去 3 ヶ月の血糖値のデータと相関する 過去 3 ヶ月の血糖値 は 調べるのが困難 しかし HbA1C は血液検査ですぐに調べられる HbA1C を調べることで 過去 3 ヶ月の血糖値が推定できる エクアドル ペルー沖の海水温 エルニーニョ現象 気候変動 気候を 予測 できる ある程度の確率で 未来が予測できる まだ見つかっていない 相関関係 がわかれば 発見 として自然科学の進歩に役立つ 教科書 P242
相互相関関数 信号 x(t) と y(t) とがどんな関係にあるのか 調べる 教科書 P242 相互相関関数の式を見てみます 二つの関数を j だけずらして 重ねて掛け合わせ 和をとる ずらして うまく重なる 時に 最大値となる関数 φ xy ( j) = 1 M M i=1 ( j = 0,1,, N M ) x(i)y(i + j)
東京の月別平均気温 http://www.data.jma.go.jp/gmd/risk/obsdl/index.php 今日はこのデータを使いますが 計算する地域 期間には指定があります 各自の該当地域 期間のデータで計算して下さい
月平均気温のデータを選択 場所を選択したら 月平均気温 を選んで下さい
期間を選んでデータを表示 期間を選んで データを表示 させて下さい 思い入れのある場所 期間がある人はそうしたデータを選んでも構いませんが 他の人と重複がないようにして下さい
Excel への貼付け Excel に貼付けて 簡単に整理しておきます 次に 太陽の南中高度を求めます http://api.knecht.jp/geocoding/q/ 東京都豊島区東池袋 帝京平成大学池袋キャンパスの緯度 :35.7328 付近 小数点以下 5 桁目は 教室が変わると値が変わります (1 歩歩くと 6 桁目 7 桁目の数字が変わります ) 太陽の南中高度 ( 度 ) = 90 - ( その場所の緯度 ) + ( 太陽の視赤緯 ) ( 国立天文台サイトより 引用 )
太陽の視赤緯 正確な数字は 以下のサイトで掲載されています http://park12.wakwak.com/~maki/de405.htm ただ 簡単のために 3 月と 9 月を 0 度として ±23.4( 地軸の傾き ) の見かけ上の運動をしているとして 三角関数で計算します A10 のセルに 帝京平成大学の緯度 (35.73) を入力し A11 のセルに =PI() という関数式を入力します A12 A16 までは 角度 ( 度 ) 角度 ( ラジアン ) 太陽の視赤緯 太陽の南中高度 と見出しを入力しておきます
関数式 B12 のセルには -60 という値を入力します C12 のセルには =B12+30 という式を入力します C12 のセルの中身を 右に 12 月までコピーします B13 のセルは =2 * $A$11 * B12 / 360 の式を B14 のセルには =sin( B13 ) の式を B15 のセルには = 23.4 * B14 の式を B16 のセルには = 90 - $A$10 + B15 の式をそれぞれ入力します
関数式の考え方 直感的にわかりやすいように 最初に 3 月を 0 度とする三角関数の起点になるように 毎月 30 度ずつ 1 年で 360 度変化する三角関数の値域を用意しました (12 行目 ) しかし 実際の三角関数の値域は ラジアン表示ですので 2π (deg) / 360 の計算で ラジアンに変換しました (13 行目 ) これで 三角関数の sin を計算しました (14 行目 ) 振幅の 23.4 度を掛けて見かけ上の太陽高度の変動 ( 視赤緯 ) を計算し (15 行目 ) 最後に 南中高度を計算しました (16 行目 ) これを 12 月までコピーします
平均気温と南中高度のグラフ 簡単のために 値だけをコピーします Excel では 一般には変化する量は行方向にとりますが 表示などの都合で 列方向に変化量をとる使い方もあります
平均気温と南中高度 形が良く似ています なんとなく 相関していると言えそうです ですが 南中高度のピークは 6 月 ( 夏至は 6 月 ) なのに 気温のピークは 8 月です これは 8 月までは 放散エネルギーよりも流入エネルギーの方が大きいため 積分されて増加している と考えられます これを 散布図で表現してみます すると 楕円になります 完全に相関している時は 散布図が直線になります 位相が 90 度ずれると 円になります
相互相関関数計算の準備 準備のために シートを変えて 行と列を入れ替えてコピーします 次に 上に 5 行挿入し 8 12 月分をコピーし下に 1 5 月分をコピーして 連続したデータにします
原因と結果 太陽の南中高度と 平均気温の関係で 気温が太陽の高度に影響を与えているのか 太陽の高度が気温に影響を与えているのかを 考えます この場合 原因は太陽の高度であり 気温が結果であると 推測します そこで 太陽の南中高度を左に 気温を右にして 準備します
相互相関関数の計算 (1) D1 のセルに -5 を入力します E1 のセルに =D1+1 を入力し これを N1 のセルまでコピーします 5 から 5 まで変化する見出しが出来ます これは 5 ヶ月から 5 ヶ月まで変化させて 位相のズレを見るためです
相互相関関数の計算 (2) D7 のセル ( 5 の列の 1 月の行 ) に 以下の式を入力します =index($b$1:$b$23, ROW() ) * index( $C$1:$C$23, ROW() + D$1) INDEX 関数は 最初のパラメータの配列から 2 番目のパラメータの番号の値を取り出す という関数です
相互相関関数の計算 (2) 意味 ROW() 関数では 現在ある行番号を取り出します 配列で 1 行目から範囲指定しているため 最初の index( $B$1:$B$23, row() ) で 現在ある行の 太陽の南中高度の値が取り出されます 2 番目の行では 現在ある行の値に D 列の 1 行目 (1 行目は絶対指定 ) の値だけ 前後にずらす計算をしています row() + D$1 の部分です 平均気温の値は この D1 のセルの値 (-5) だけずらして取り出します ですので 平均気温の項は index( $C$1:$C$23, row() + D$1 ) になります
数表を埋める 式は 絶対参照と相対参照を組み合わせて作っていますので 一気にコピーして 面を埋めて下さい
相互相関の計算 数表を計算したら D24 に =D1 の式を書き D25 に =average(d7:d18) の式を書きます これを N 列までコピーします
空白行の意味 セルの上下に空白がある意味を考えてみて下さい 過去 5 ヶ月分とか 未来の 5 ヶ月分との関係を調べます 比べる相手が上下に 5 ヶ月分ずつ幅広くなっているので 時間で言えば前後 5 ヶ月分は 相手のデータ は必要でも 自分自身のデータ は使わない領域があります その領域は 計算式だけ入れても答えは出てきません
相互相関関数のグラフと読み方 24 行目を軸ラベルとして 25 行目の数値だけをグラフ化します 位相 (5 ヶ月前 ) (5 ヶ月後 ) まで ずらして計算した結果 相関が最大となるまでに 1 2 ヶ月の位相のズレがあることがわかります
位相補償 前処理で こうした位相のズレを補償すると 二つの変化が 一次関数 の形に近くなります ( 直線関係 ) 実際に 平均気温と 南中高度のデータを 2 ヶ月ずらして散布図を作成すると かなり 右肩上がりの直線に近づきます 良い相関関係 があるデータでも 時系列データの場合 位相をきちんと補償してあげないと 相関係数 は悪くなる場合があります 月 単位ではなく 日 単位で計算して位相保証すると 相関はほぼ 1 になります
相関係数の計算 相関係数は 以下の式で計算できます n n (x i x)(y i y) i=1 (x x) 2 (y i y) 2 i=1 n i=1
相関係数の考え方 x や y が 平均よりどの程度大きいか 小さいか ペアになる 1 項目ごとに両辺がそれぞれ x y の長方形を作ってみます 傾向が一致していれば x が大きいときは y も大きく 合計はプラスになります マイナスとマイナスを掛けてもプラスになるので 傾向が一緒だとプラス側だけに長方形が来ます 傾向が真逆だとすると x が + の時は y がー x がーの時は y が + となって 合計はマイナス側で絶対値は大きくなります 特定の傾向がない場合は 相殺して 0 に近づきます この長方形の面積を x や y の大きさで揃えると ( 分母の式で割ると ) 最大が 1 最小が 1 になります この分母は標準偏差が使われます
相関係数の式変形 この式のままだと 一回平均を計算しないと相関係数が計算できません そこで x = 1 n n x i, y = 1 i=1 n n i=1 という関係を利用して 式変形します ( 板書参照 ) 最終的には 以下の式になります y i r = n i=1 n i=1 x 2 i 1 $ & n % x i y i n i=1 x i 1 n ' ) ( 2 n n x i y i i=1 i=1 n i=1 y 2 i 1 $ & n % n i=1 y i ' ) ( 2
Excel での相関係数計算の準備 Excelの機能では 簡単に相関係数を求めることもできますが あえて 式の意味から考えて欲しいため バラバラに計算します n 右の式を用いて Σ を置き換えると 相関係数は下記のようになります r = S XY 1 n S X S Y S XX 1 n S 2 X S YY 1 n S 2 Y S X = S Y = S XX = S YY = S XY = i=1 n i=1 n i=1 x i y i n i=1 n i=1 x i x i y i y i x i y i
Excel での相関関数 2 ヶ月分だけ位相をずらします そこで X Y XX, YY, XY を各行ごとに計算します ( 具体的な手順は 先週のスライドを参照して下さい ) それらの値から r を求めると r=0.9592 が得られます
相関係数の読み方 一般的に 例えば医療関係のデータを読む際には 0.7 1: 強い相関がある 0.3 0.7: 相関を伺わせる傾向が見られる -0.3 0.3: 無相関である -0.7-0.3: 負の相関の傾向が見られる -1-0.7: 負の相関がある などという言い方をします 0.3 から 0.7 の間は 非常に扱いが難しいことが多い 但し 相関係数は X と Y を入れ替えても同じ結果になるので どちらが原因で どちらが結果かを言うことはできません
とっても楽な計算方法は 実は =CORREL(B2:B13,C2:C13) と入力すると 相関係数は計算できます ただ 元となる考え方や 感覚的に 何を計算しているのかを感じ取ってから 数式を理解し 最終的には =correl() を使うようにして下さい そうしないと 位相差のあるデータ をそのまま分析して 目立った相関は見られない などという論文を書いたりしてしまいます ( 後で 結構恥ずかしいのではないかと思います )
自己相関関数 今度は 関数の波形を自分自身と相関させてみます 第 10 回で使った 心電図のデータについて 相互相関の代わりに 自分自身の値との相関を出してみます ( 私は 10 刻みで位相を変えました そうしないと データが大きいので )
計算を手動に切り替え Excel で サイズの大きいデータの自己相関を計算させると 何度も 再計算 が実行されて 非常に 重く なります 何か入力するたびに 全部を再計算するために 反応が鈍くなります そこで 一度数式を入力して計算させたら 自動 計算をやめて 手動 再計算に切り換えます または ( 裏技ですが ) 計算済みのセルの値を 全部コピーし 値を貼付け して 計算しないようにしてしまいます こうした 裏技 を使うときは 自分のやったことをしっかり覚えておくように
自己相関の意味 自己相関の結果を見ると -150-170 の付近で値が大きくなっています これは 150 サンプル 170 サンプル (0.75 秒 0.85 秒 ) で データが反復していることを示唆しています これの山は 心電図の周期を意味しています が 計算する区間が広いと 大まかな値 しか出てきません ( 瞬時心拍数の算出には不向き )
自己相関と周期性 周期性のある信号では 自己相関を計算すると その周期に相当する 遅れ の部分で 相関関数がピークを持ちます ということは 自己相関を計算すると 信号の周期性が計算できるのではないか? と考えられます これが ウィーナー = ヒンチンの定理 (Wiener Khinchin theorem) と言われ 自己相関関数のフーリエ変換は 信号のパワースペクトルになる という関係式です この科目の守備範囲外ですが フーリエ変換との関係は有名ですので 覚えておいて下さい
今日の提出課題 太陽の南中高度と平均気温のデータを用いて 相互相関関数のグラフを作成して下さい (2 点 ) 5 ヶ月から +5 ヶ月まで 平均気温のデータをシフトさせて 南中高度と平均気温の相関係数を計算し 何カ月ずらした時に相関係数がいくつだったか 表にしてまとめなさい (3 点 ) ヒント :INDEX 関数をうまく使うと 一気に作成できます CORREL 関数を使っても構いません この結果をグラフに表示して下さい (1 点 ) 以上のレポートを 今日の 出席点 とします 出席そのものは カードでチェックします
発展課題 ( ボーナス課題 ) 気象庁のデータでは 週単位 の区切りのデータもダウンロードできます できれば 10 年分程度を平均したデータを使って下さい 太陽の南中高度のデータも 週単位 で計算し 平均気温と南中高度のデータの 位相差 を求めて下さい 夏至 から 何週間くらい遅れて 気温のピーク が来るか 体感している時期が結果として得られているか 検証して下さい 簡単のため 1 年を 52 週間 (364 日 ) として構いません ( ボーナス点 : 最大 4 点 )
次回の予告 第 5 回 方式 性能 誤差 位法と零位法 感度 分解能 測定範囲 スパン ダイナミックレンジ 直線性 確度 応答時間 分散 標準偏差 標準誤差 有効数字について学習する 計測論 の部分を学習します
もう一つの提出課題 欠席課題 今日の講義ノートの 欠席課題 ですが 授業課題ではなく こちらの課題で出しても構いません 相関があるということは何を意味するかを最初に説明し 相関係数を求める計算式で なぜ相関が表せるか 数式に基づいて説明を書いて下さい この課題についても 6 点満点で採点します 授業課題のレポートで提出しても構いませんし こちらの課題で出しても構いません
欠席した人は 今日の授業課題を提出して下さい 提出があれば 出席に切り換えて レポートも採点評価します