情報処理実習2

今回はデータ処理の一つの要である最小二乗法について学ぶ．

2つの変数間の関係が理論やモデル等で判明しているものを実験などで数値データを計測し，その係数を求めたい場合がある．
（例えば比例関係となるばね係数や，1次遅れ系の時定数の決定など）

例えば，比例関係であれば最小で2点の計測データから直線を1つ決めることができるが，計測データなどには「測定誤差」が含まれるため，より正確な係数を求めたい場合は必要なデータ数よりも多めに計測し，それらを全て用いて最もよくあてはまる直線を決めると良いとされている．
この際，得られたデータ群をもとに客観的な方法で最もよく当てはまる関数の係数を決める方法の代表例が，最小二乗法である．
最小二乗法で決定できる関数の形は様々であり，工学でよく用いられる関数は１次関数の他に，２次関数，高次の多項式，指数関数，周期関数など，様々な場面で用いられる．

補間法が「制御点を滑らかな曲線で繋ぐ」ことで，制御点以外の関数値を求めるための手法であるのに対し，最小二乗法は関数形をあらかじめ与えて，計測などにより得られた点群に最もよく当てはまる関数の係数を決定することを目的としており，その目的が大きく異なる．

最小二乗法の考え方

いま，比例関係にある計測データとして $N$ 個の点群 ($x_i , y_i$) が得られたとする．
これらをもとに，もっともよく当てはまる関数（この例では直線）の係数を決定する方法が，最小二乗法である．

ここで，求めたい直線の式を $y=ax+b$ とする．最小二乗法は，この $a, b$ の値を決定することが目的である．

直線と計測点との差（図中の青色矢印の部分）は，当てはめの度合いを示しており，残差と呼ばれる．
残差の合計が小さいほど，直線が点群の傾向をよく表していると言える．
各点における残差の値は正にも負にもなるので，そのまま足すと正負でキャンセルされてしまうためよろしくない．そこで，残差の二乗の和を最小化することを考える．

残差の二乗和を $e$ とすると，$e$ は直線の傾きと切片 $(a,b)$ の関数となる．
以下の式で表される．

\begin{align} e(a,b) = \sum_{i=0}^{N-1}\left( y_i - (ax_i+b)\right)^2 \end{align}

$e$ は $ a, b$ の2次関数となるので，これが最小となる必要条件は，

\begin{align} \dfrac{\partial e}{\partial a} &= \sum_{i} 2\left( y_i - (ax_i+b)\right)(-x_i)=0 \notag\\ \dfrac{\partial e}{\partial b} &= \sum_{i} 2\left( y_i - (ax_i+b)\right)(-1)=0 \end{align}

\begin{align} \sum x_iy_i - a\sum x_i^2 - b\sum x_i=0 \notag\\ \sum y_i -a\sum x_i -b\sum 1=0 \end{align}

となる．これを$a,b$ に関する連立一次方程式として，$a, b$ について解くと

\begin{align} a &= \dfrac{\sum x_i\sum y_i - N\sum x_iy_i}{\left(\sum x_i\right)^2-N\sum x_i^2} \notag\\ b &= \dfrac{\sum x_iy_i\sum x_i - \sum x_i^2\sum y_i}{\left(\sum x_i\right)^2-N\sum x_i^2}\\ \end{align}

となる．
一見複雑に見えるが，それぞれ右辺は既知の点 $(x_i, y_i)$ の算術演算のみで表されるから，簡単な代数計算により係数 $a, b$ が求まる．

ここでは，簡単な直線の例で示したが，2次関数や指数関数の最小二乗法も，残差を最小にするという考え方は同じである．

課題1

実験で，あるばねの変位と荷重に関するデータを得たとする．
以下のソースコード内の配列 x[], y[] はそれぞれバネの変位と荷重のデータである．
荷重と変位の間には比例関係があるとして，この比例係数（ばね定数）を求めよ．

まず，ソースコード中のコメントを参考にして，式(4)右辺に出現する各値（$x,y$ の和，和の二乗，二乗の和，内積など）を計算せよ．関数にしておくとよい．
次に，最小二乗法による一次関数の係数（直線の傾き $a$ と切片 $b$ ）を求めよ．
計算結果の値を画面に表示するプログラムを作成せよ．

参考ソースコード：

#include <iostream>
using namespace std;

double sum(double x[], int n)
/* 配列の合計 Σ を計算する関数 */
{
    ...
}

double sum2(double x[], int n)
/* 配列の2乗の合計 Σ を計算する関数 */
{
    ...
}

double inner_product(double x[], double y[], int n)
/* 二つの配列 x,y の内積を計算する関数 */
{
    ...
}

int main(void)
{
    const int N = 10;

    // あるばねにおける変位 x[mm] に対する荷重 y[N]
    // 値を少し変化させてみて，各自で結果の変化を確認せよ．
    // 動的配列を使用してもよい．
    double x[N]={0, 1.1, 2.1, 3.0, 3.9, 5.0,  5.9,  7.1,  7.8, 8.9};
    double y[N]={0, 2.4, 4.2, 5.7, 8.2, 9.2, 11.9, 14.0, 16.2, 17.4};

    double a, b;    // 一次関数の係数

    //  以下にコードを記述

    //  結果の出力
    cout << "a = " << a << ", b = "  << b << endl;

    return 0;
}

実行例：(値が正しいとは限らない)
a = 0.977565
b = 0.063278

ソースコード 153r000000-??-1.cpp を提出せよ．

課題2

上記の計算コードを利用して，計算により以下のようなグラフを描け．
（この図では，直線と点群の最適なフィッティングができていない．）

グラフの作成には，Excelやグラフ作図ソフトを使用しても良い．または，gnuplotを使用することもできる．

参考までに，gnuplotで作図するテンプレート springrate.cpp を利用してもよい．
Gnuplotによる作図方法についての詳細は，本ページ末尾を参照．

ソースコード 153r000000-??-2.cpp およびグラフを画像形式(jpg, png, pdfなど)で保存したもの 153r000000-??-2.jpg, png, pdf を提出せよ．

グラフの作成法

プログラムによる計算やデータ処理の結果は，数値の羅列であるため，視覚的に表示することはとても重要となる．
それだけでなく，実験レポートや学術論文では，グラフの書式について一定のルールがある．

グラフの作成にあたっては，表計算ソフトではなく専用のツールを用いた方が容易であり，作業効率も良い．

gnuplotによるグラフ作成

gnuplot(グニュープロット，グヌー？　ヌー？)は，昔から利用されているグラフを描くためのソフトウエアであり，無料で利用できる．
情報処理教室のPCには既にインストールされている．

このように，簡単なコマンド入力で様々なグラフを描画することができる．
もちろん，関数形だけでなく，プログラムで出力した数値データのグラフ化もできる．

方法1：Cプログラムから数値データを一旦ファイルに出力し，gnuplotでデータファイルを読み込む

自作プログラムで何かの計算を行い，データファイルとして一旦 .txt や .csv ファイルに配列データなどを保存したのち，gnuplot を起動してコマンドからファイルを読み込ませる方法である．
この方法は，データファイルを保存しておくことができるため，データ処理プログラムとグラフ作成（可視化）過程を分けることができるメリットがある．

例えば，データファイル data.csv を出力したとする．これは，さきのバネの試験データとする．これをgnuplotで読み込んでグラフ化する．
（先頭の # の行は gnuplot ではコメントとして無視される．）

まず，作業フォルダに，data.csv をダウンロードし，[File][Change Directory]で.csvファイルのあるフォルダを指定する．
次に，gnuplot の入力プロンプトで，以下のように打ち込んでみよう．

たったこれだけで，.csv ファイル内に書かれた x,y 値の散布図グラフを書くことができた．

さらに，gnuplotに与えるコマンド（この2行）をあらかじめ .plt ファイル（例えば，ex1.plt）に保存しておくと，

とするだけで，プログラムのように呼び出すことができる．

方法2：Cプログラム中から，gnuplotを直接呼び出して描画する

この方法はCプログラム中から直接 gnuplot を起動し，作図コマンドを直接送信することでグラフをリアルタイムに表示する方法である．
計算中のデータをダイレクトに可視化できるため，時間のかかる計算処理の途中結果を逐一チェックしたり，また，データ保存が必要なくグラフだけあれば良い場合に効率が良く，便利である．

Windows OS やUNIX には，あるプログラムから，外部の別のプログラムを呼び出し（実行し），データの受け渡しをする便利な「パイプ」と呼ばれる機能がある．
これを使えば，自分の作成したプログラムから gnuplot を実行し，計算データと各種設定を文字列として gnuplot に送信し，瞬時にグラフを描くことができる．

具体的には，Cプログラムから popen() 関数を使って gnuplot 実行ファイルを呼び出し（=実行），グラフ作画に必要なコマンドとデータをパイプ経由で送信する．

推奨環境以外でのパイプ処理上記ソースコードについて，本授業の推奨開発環境以外(bcc32など)では，以下のような変更が必要となる場合がある．

popen() を _popen() に，(_ を追加)
pclose() を _pclose() に，(_ を追加)
gnuplotのパス指定，例えば，const char gp_cmd[] = "/Applications/gnuplot.app/bin/gnuplot -persist";など

これ以外にも，PC環境・設定によって追加の操作が必要となる可能性がある．
自己解決できない場合は，まず情報処理教室のPCを使用して動作を確認すること．

データのグラフ表現：最小二乗法

最小二乗法の考え方

課題1

課題2

グラフの作成法

gnuplotによるグラフ作成

方法1：Cプログラムから数値データを一旦ファイルに出力し，gnuplotでデータファイルを読み込む

方法2：Cプログラム中から，gnuplotを直接呼び出して描画する