計測工学 I 第 3 回 Excel による回帰式の計算
今日の内容 第 3 回 Excel による回帰式の計算 シラバスより 第 3 回 回帰式の計算 Excel を用いて測定データから最小二乗法によって 回帰式の計算を行う この計算方法を学び 実際のデータに適用して回帰直線をグラフ化する 最小二乗法によって 計測データが満たしている関数式を推定する方法を学びます
回帰式とは何か? 教科書 P255 計測データには 誤差 がつきものです 誤差を 測定時 にどう減らすか センサの装着方法などを工夫する 差動入力回路を用いる 第 6 回 測定されたデータに残っている 誤差 をどう処理するか 統計的に扱う 標準偏差や 平均を調べる 後期第 13 回 関数形を推定する データの持っている特性を調べる 推定した関数に実際のデータを適応させ 係数を求める 今日のテーマ
今日行うこと Excel 上で 一次関数を作る ランダムな雑音を乗せる 雑音が乗った一次関数のデータから元の一次関数を推定する 二人一組 または 3 人一組で演習を行います それぞれが 雑音入り のデータを作成し それをもう一人に提示します データ点数は 21 点 : X { -10, -9, -8, -1, 0, 1, 2, 9, 10} 相手に画面を見せて 入力してもらって下さい ( または メール添付 ) データを受け取った側は そのデータから 元になる一次関数を推定する
最初の練習 (Step 1: 元の一次関数 ) Excel を立ち上げて下さい y=ax+b の関数形を 適当に決めます ここでは例題として y=2x-3 を使ってみます Sheet1 の A1 に X B1 に Y と見出しを入力します A2 のセルに -10 を入力し A3 のセルに =A2+1 と式を入力します A3 のセルを オートフィルで A22 まで下に引っ張ります X {-10, -9,.., -1, 0, 1,.., 10 } の出来上がり!
練習 (Step 1: 元の一次関数 2) B2 のセルに 関数式 =2*A2-3 を入力します B2 のセルをオートフィルで B22 まで下に引っ張ります 前ページでも使いましたが オートフィルは セルの右下の の上にポインタを重ねると十字カーソルに変化するので その時に右や下に移動させると実行できます 元の関数値の出来上がり
出来た一次関数をグラフで確認する A1:B22のセルを選択します 挿入 グラフ 散布図 を選びます 散布図は 横軸はX 縦軸はYというグラフです 平滑線を選ぶと良いでしょう
雑音を載せる ( 乱数生成 ) C1 のセルに 乱数 と書きます C2 のセルに =rand() と入力します C2 のセルの式を C22 までオートフィルでコピーします
乱数の範囲を変える Excel の rand() 関数は 0 rand()<1 の範囲の乱数を合成します 値域 [0 1) を [-1, 1] に変換したいので D 列で =2*C2-1 という一次式を入力します そうすると D 列には 1 から 1 の間の乱数が作られます ちなみに a x<b,a<x b をみたすの集合をそれぞれ [a,b),(a,b] と表わし 左閉右開区間 左開右閉区間という (http://www.suriken.com/knowledge/glossary/half-opening-section.html) 数式当てゲーム をやります 相手に意地悪するために 乱数を大きく [-2,2) や [-3,3) の乱数にするには どんな式にしますか? ただし 乱数の中心は 必ず 0 になるようにして下さい 自分で自分に意地悪するのでも構いません
雑音レベルを変える 一次関数は写像です [0, 1) の値の範囲を例えば [-3, 3) に変えるには y = ax + b で x=0 の時 y=-3, x=1 の時 y=3 となるように a と b を求めれば良いことになります ( 中学校 1 年生あたりで習ったはず ) -3 = a 0 + b, 3 = a 1 + b を連立させます そうすると b = -3, a = 6 が得られます つまり y=6x-3 D 列のセルに =6*C2-3 と入力すれば [-3, 3) で変化する乱数が得られます
乱数の乗った関数値 グラフを簡単に作れるように E 列に X 列の値をもう一度表示します E1 のセルに X を E2 のセルに =A2 を入力します F 列に 乱数の乗った関数値を用意します 但し 桁数が多いので小数点以下 3 桁に丸めます F2 のセルに =round(b2+d2,3) という式を入力します これを F22 までコピーし 散布図にしてみます ( コピーの方法や散布図の作り方は説明済み )
誤差と測定値 今回作成したデータ列は 一次関数です そこに誤差 ε が乗って来ていますので 式で表すと y=ax + b + ε で表すことができます 誤差 ε でかく乱されていますが 測定値 x, y の間には a, b という 測定対象の性質 が隠れていて これを探し出すことが 測定 の目的です その 真の姿 をかく乱する ε を ε = y (ax + b ) と置いて この影響を取り除ことを考えます
係数 a, b の推定 もし 係数 a, b がわかっているなら ε = y (ax + b ) は 計算で求められますが わからないので推定します X と Y は 実測値 がわかっていますから a や b を適当に当てはめてみると ε が求まります 適当に a や b を代入して実測値から ε を計算して 誤差が最小となる a, b の組合せが 真の係数 である と考えます ですが 適当に a, b を当てはめてもうまく誤差が最小になるとは限りません 教科書 P256
最小二乗法の考え方 ε = y (ax + b ) の両辺を二乗します 誤差は ± に振れますが 二乗すると 0 を中心として左右に増加する関数になります これを 全ての測定点について加算します 仮に この式を Ε と置きます Ε = ε 2 = {y (ax + b)} 2 ここで a を変化させた時 a が真の値と一致した時にこの式は最小となります 同様に b を変化させた時に b が真の値と一致した時に この式は最小となります
a と b による偏微分 a を変化させていった時 a が真の値の時 ( これを a t とします ) E は最小になります 同様に b の真の値を b t として b が b t の時に E は最小になります Ε = ε 2 = {y (ax + b)} 2 つまり E を a や b で偏微分すると 変曲点である真の値で傾きが 0 となりますから 以下の式が成り立ちます a Ε(a ) = 0, t b Ε(b ) = 0 t
右辺の展開 右辺について x や y は実測値で 既知 です a や b がこれから求める係数です ですので a や b に注目して 式を書き直します {y (ax + b)} 2 = {y 2 2y(ax + b) + (ax + b) 2 } b 2 の項は x や y を含まないので 定数 として扱える つまり Σ で加算した回数 N 倍すれば良い = (y 2 2axy 2by + a 2 x 2 + 2abx + b 2 ) = y 2 2a xy 2b y + a 2 x 2 + 2ab x + Nb 2
右辺の展開 (2: 定数の Σ) 定数の Σ b 2 は x や y に対して変化しない値 つまり定数です 定数を Σ する場合 (b 2 + b 2 + b 2 + + b 2 ) と Σ で加算する回数分だけ加算するため Σ で加算する個数を乗じれば良い ということになります 今回の Excel のデータの例では { -10, -9, -1,, 0, 1,.. 10} と 10 から 10 までの範囲で X を変化させていますから N=21 になっています b 2 = Nb 2
右辺の展開 (3: 測定値の扱い ) X や Y は 実際に測定した値になります つまり 全て既知です これに対して a や b を変化させて関数 Ε を調べるので X や Y の Σ は 全て既知の測定値を用います そこで 右のように置いて 式を書き直します また この式の変数は a や b ですので S は定数として扱います S XX = x 2 S YY = y 2 S XY = S X = S Y = x y xy y 2 2a xy 2b y + a 2 x 2 + 2ab x + Nb 2 = S YY 2S XY a 2S Y b + S XX a 2 + 2S X ab + Nb 2
Ε の a による偏微分 変形した式を a により偏微分します b も変数ですが a で偏微分するときは b は定数として扱えます a が真の値となった時に この式の値が 0 になります Ε = S YY 2S XY a 2S Y b + S XX a 2 + 2S X ab + Nb 2 a による偏微分では a を含まない項はすべて 定数 として扱える 定数の微分は 0 なので a を含む項だけ考える a Ε = 2S XY + 2S XX a + 2S X b a や b が真の値の時に この式が 0 となる S XY + S XX a + S X b = 0
Ε の b による偏微分 同様に b により偏微分します a も変数ですが b で偏微分するときは a は定数として扱えます b が真の値となった時に この式の値が 0 になります Ε = S YY 2S XY a 2S Y b + S XX a 2 + 2S X ab + Nb 2 b による偏微分では b を含まない項はすべて 定数 として扱える 定数の微分は 0 なので b を含む項だけ考える b Ε = 2S Y + 2S X a + 2Nb a や b が真の値の時に この式が 0 となる S Y + S X a + Nb = 0
求める a, b の連立方程式 ここで 二つ変数 a と b について 二つの式が得られました S XX a + S X b S XY = 0 S X a + Nb S Y = 0 この連立方程式を解くと 解は以下のようになります a = NS XY S X S Y NS XX S X 2 b = S XX S Y S X S XY NS XX S X 2
Excel での Σ 計算 作成した X と Y の系列のグラフを 同じテーブルなどの友人に メールで送るか 見せて手入力してもらって下さい (USB でのコピーでも OK です ) まず 元の関数式を隠した形で X と Y の値のデータだけ 新しいファイルにコピー します Excel で 新しいファイルを作成し そこに 数値 だけをコピーし 作成したその新しいファイルを友人に送ります
出題の手順 Step 1: 新しいファイル を作成する Step 2: 出題用に B 列 (Y) の値の式を全部書き換える Step 3: 出題用に ノイズレベルを変える Step 4: 作成した X と Y のデータ (E1:F22) を新しいファイルにコピーする この時 値だけを貼付け する Step 5: 名前を付けて保存 し メールで友人に送る
測定データの確認 友人が作成した ノイズに埋もれた一次関数の数表を あなたが 測定したデータ と見なします この数列から 測定データを推定します ここで 以下の五つのパラメータを計算します Σ を求める前に まず 各 x, y ごとに x 2, y 2, xy を計算します C 列 D 列 E 列に計算結果を入れることにします S XX = x 2 S YY = y 2 S XY = S X = S Y = x y xy
各列の計算式 C1 は X*X とします D1 は Y*Y とします E1 は X*Y とします C2 のセルには =A2*A2 の式を書きます X*X, X 2 です D2 のセルには =B2*B2 の式を書きます Y 2 です E2 のセルには =A2*B2 の式を書きます X*Y です C2:E2 のセルに 式を入力したら 三つのセルを選んで 一気に C22:E22 まで オートフィルでコピーします 実データから N を求めるため F1 のセルに =count(a2:a22) の式を入力します N=21 になるはずです
Σ の計算 各行の値が求まったら 合計 (Σ) を計算します A 列で縦方向に 全部の要素を選びます 選んだら オートサムで合計をクリックして下さい または =sum(a2:a22) の式を入力して下さい 同じ作業を繰り返すか A23 のセルをオートフィルで E23 までコピーします
a と b の推定値を求める A23 には S X が B23 には S Y が C23 には S XX が D23 には S YY が E23 には S XY が 計算されています また 今回の例では N=21 が F1 に計算されています さすがに この求め方は 各自で工夫して下さい 私の例では 雑音レベルを 3 にして a=1.06, b=-3.29 になりました a = NS XY S X S Y NS XX S X 2 b = S XXS Y S X S XY NS XX S X 2
回帰式との同時表示 EXCELの散布図では 自動計算もできますが 測定値の列に並べて計算値を表示させ 比較することができます 3 列で散布図を作ってみて下さい 計算で求めた a の値 b の値は $F$2 など セルの絶対アドレスで指定します 右の例では I 列に =$F$3*G2+$F$4 と書いています
回帰式を求める意味 計測したデータでは 雑音が含まれているため 元の関数形を そのままでは求めることが困難です 元の関数形を調べることにより 測定を通じて 対象物の 特性 などを知ることができます 今回は 一次関数で解きましたが 教科書 P258では2 次関数の例があります さらに exponential でも 対数を取るなどして近似ができます ( が 今回は扱いません )
今日のレポート (1) 自分で 最小二乗法の係数を求める計算をして欲しいと思います ( 各自が考えて下さい ) レポート課題に含めてもいいのですが きっと 式を丸写しするだけで全然頭を働かせずにレポートを書くと思うので そんな 採点が面倒なだけの課題は出したくありません 友人から 出題 された ノイズだらけの一次関数 のデータを処理し a と b の値を推定して下さい その計算結果に 友人の作った 正解 を書き添えて [ 学籍番号 -03] のファイル名で メール添付で提出して下さい 今日の 出席 は カードで見ます ここまでで 6 点
今日のレポート (2) 係数を求めたシートで 問題の X と Y の列に並べて 計算で求めた 回帰式 によるデータ系列を計算し X に対する Y( 測定値 ) と Y( 回帰式 ) のグラフを表示して 作成して下さい このグラフで 4 点を加点します 今日のレポートの配点は 10 点です
今日のレポート ( ボーナス課題 ) 教科書を参考にして 友人が作成したデータに重畳させた雑音のレベルがどの程度だったか 測定値から推定する方法を 式を省略せずに求めて 報告して下さい 数式ベースで示すことが基本です 考え方や結果と合わせて 内容が正しいレポートには最大で 5 点まで加点します
次回の予告 第 4 回 Excel による相関係数の計算 シラバスより以下の部分を扱います 第 4 回 Excel による相関係数の計算 相関の考え方を導入する 相関が高いということ 負の相関があるという概念を学び 時系列データの位相差 相互相関関数 相関係数 自己相関関数について学ぶ 二つのデータが伴って変わる時に どんな関係が成り立つかを求めるのが相関係数です 次回も今日のシートを使います 削除しないでおいて下さい 削除した場合は 今日の資料に沿って また作り直して下さい