equal_l2’s blog

※記載されている内容の正確性は保証しませんが、間違いを指摘していただければ直します。

最小二乗法(7) 項の欠けた線形近似

物理・数学

前:最小二乗法(6) 線形近似における係数の推定誤差 - equal_l2’s blog
次:最小二乗法(8) 重み付き多項式フィッティング - equal_l2’s blog

表題がどうも上手くつけられなかったので、シチュエーションを説明しよう。

今までやってきた\ d\ 次の多項式近似で使っていたモデル関数は、0次から\ d\ 次までの項が必ず1個づつあって、欠けることはなかった。

しかし、時には「0次の項はいらなくて、1次の項だけ欲しい」などというシチュエーションもあるかもしれない。
この場合、単に1次の多項式近似を行って0・1次の項を推定したうえで0次の項を無視する、ということはできない。0次の項なしには残差の二乗和を最小にはできないからだ。
(推定の仕方からして当然である)

こういう時はどうするかというと、偏微分で書いた正規方程式まで戻るしかない。

\def\vb{\vec{\beta}}\def\fxb{f(x_i,\vb)} \displaystyle \sum_{i=1}^{n}\frac{\partial \fxb}{\partial \beta_j} y_i = \sum_{i=1}^{n}\frac{\partial \fxb}{\partial \beta_j} \fxb
誤差なども、すべて一から計算し直しである。

モデル関数f(x)=\beta_1 xについての線形近似

結果を先に申しあげると、係数は
\displaystyle \beta_1 = \frac{Y_1}{X_2}
その誤差は、
\displaystyle\begin{eqnarray} \sigma(\beta_1) &=& \sigma_y \sqrt{ \frac{1}{X_0 X_2} } \end{eqnarray}
となる。
ただし、\sigma_y=\displaystyle{\sqrt{\frac{\sum_{i=1}^{n}\bigl(y_i- \beta_1 x_i\bigr)^2}{n-1}}}

続きを読む

横に動ける二重振り子(1) ラグランジアンを立てる

解析力学

友達からこんな物理モデルの運動をアニメーションできないかと言われた。
\def\to{{\theta}_1} \def\tt{{\theta}_2} \def\lag{{\cal L}}
f:id:equal_l2:20150420212540p:plain
振り子の下にもう一つ振り子をつなげた、いわゆる二重振り子というやつである。
ちょっと普通と違うのは、この振子は一番上がx方向に動くことである。
(x軸という棒に輪っかを通して、そこに二重振り子を結びつけたと考えればいいだろう)

何はともあれ運動方程式を立てないとどうにもならないので、立ててみる。

二重振り子の問題はラグランジアンで解くと相場が決まっているので、ラグランジアンを考える。
一般化座標として、一番上の点のx座標であるx、上の振子と鉛直方向がなす角\to,下の振子と鉛直方向がなす角\ttの計3つ(図参照)を用いる。

結論を先に書いてしまうと、
\begin{eqnarray} \lag &=& \frac{M}{2} \left( \dot{x}^2 + 2 l_1 \left( \dot{\to} \dot{x} - g \right) \cos{\to} + l_1^2 \dot{\to}^2 \right)+ \\ && \ \ \ \ \ \frac{m_2}{2} \left( l_2^2 \dot{\tt}^2 + 2 l_2 \left( \dot{\tt}\dot{x} - g \right) \cos{\tt} + 2 l_1 l_2 \dot{\to}\dot{\tt} \cos{(\to - \tt)} \right) \end{eqnarray}
この後はつらく苦しい導出過程。

続きを読む

参照型の引数について

C/C++ プログラミング

関数の引数を参照型でとることには、ポインタを使うのと比べていくつかの利点がある。

1.呼び出す際にアドレスを渡さなくていい

変数の交換を行う関数でaとbを交換することを考える。
その関数をswapとする。

ポインタで実装したswap関数は、次のようにつかう。

swap(&a,&b);

一方で、参照で実装したswap関数は、次のようにつかう。

swap(a,b);

違いは明確で、ポインタ実装ではアドレスを渡すのに対し、参照実装では変数(オブジェクト)を渡している。
ひと手間減るわけだ。

2.書き方を変えないで済む

int型に対して、ポインタで実装したswap関数はこんなようなのだった。

void swap(int* a, int* b){
 int c = *a;
 *a = *b;
 *b = c;
}

一方で、参照で実装すると

void swap(int& a, int& b){
 int c = a;
 a = b;
 b = c;
}

参照実装のほうは、関数の中身をそのままmain関数などに貼り付けても問題なく動作する。
ポインタ実装のほうはそうはいかない。

これらの利点は、要するに「変数に直接アクセスする」ことを楽に行える、とか意識せずに行える、ということである。
逆に、変数が変更される可能性があることに気付きづらくなるので、企業のコーディングスタイルではconst属性の参照(参照先の値を変更できない)以外を引数の型にすることは禁止されていることもある。

参照を学ぶ

C/C++ プログラミング

ポインタとよく一緒に語られる存在として、C++の参照がある。

ただC++の参照というのはオブジェクトに別の名前を付けることらしく、実装にポインタを必ずしも使う必要はないらしい。

参照型変数

参照型の変数は、次のように宣言する。

参照先の型 &変数名 = 参照先オブジェクト;

たとえば、int型の変数aへの参照であるraを宣言するには…

int& ra = a;

とする。

参照型の性質として、次のことがあげられる。

1.宣言時に必ず初期化する必要がある。

宣言する際に、必ず参照先を指定しなければならない。
つまり、参照型の変数は宣言のみを行うことはできない。

2.参照先は変更できない

ポインタと違い、示す先、もしくは参照先を変更することはできない。

3.参照先の変数の代わりとして使える

参照型変数は、参照先オブジェクトと全く同じように扱える。
(ポインタのように間接演算子でアクセスする必要はない
もっとも、参照型変数はWindowsのショートカットとかUnix系のシンボリックリンクのようなものなので当然といえば当然なのだが。

参照型引数(参照渡し)

参照型は、ポインタ同様に仮引数の型としても使える。
しかも、ポインタのようにアドレス演算子を使って引数を渡さなくてもいいのが利点となる。
これを使えばswapの素朴な実装が動くようになる。

#include <cstdio>
void swap(int& a, int& b){
 int c = a;
 a = b;
 b = c;
}

int main(){
 int a=10, b=20;
 printf("before:%d %d \n",a,b);
 swap(a,b);
 printf("after :%d %d \n",a,b);
}

(C++なんだからほんとはテンプレートで型汎用的に書いたほうがいいんだろうけど)

ただし参照型引数では、引数が関数によって変更されるのかどうかが分かりづらいという欠点もある。
実際 google のC++スタイルガイドでは、引数を変更する際にはポインタ型引数として取るようにせよとの記述がある。

実はあんまり使い道ない?

ポインタを学ぶ(4) おまけ

C/C++ プログラミング

あんまりポピュラーではないけど、ちょこちょこ出てくるテクニック。
（ポインタともあんまり関係ないけど）

条件分岐文の節約

0以上の整数入力に対して、1個づつ異なった値を、結果として返す処理を書いてみる。

以下のコードで共通の変数として、整数入力を格納する変数をint in;とする。

とんでもなく愚直に書くなら、if文を使って…

int in; //入力を受け取る変数
int data0 = ? ,data1 = ? ,data2 = ? , ...
...
if(in == 0){out = data0;}
else if(in == 1){out = data1;}
else if(in == 2){out = data2;}
...

Cを学ぶ上で割と最初に知るべきこととして、同じようなものを入れる変数は配列にまとめたほうが良い、ということがある。
この原則をdata0, data1, data2, ...に適用して、data[]にまとめる。
ただし、配列のインデックスと入力は対応するようにする。
たとえば、in == 0に対してはdata[0]が対応するようにする。
こうすると、次のように書き換えられる。

int in; //入力を受け取る変数
int data[] = { ... };
...
if( in >= 0 && in < sizeof(data)/sizeof(data[0]) ){
 out = data[in];
}
...

配列を扱うので、添え字に問題がないか(正の数か、要素数を超えていないか)をチェックする必要はあるものの、処理そのものは少なからずスマートになった。

ちなみに、sizeof(data)/sizeof(data[0])は、配列の要素数を求める非常にポピュラーな記法である。
標準Cではこの方法以外に配列の要素数を求める方法はない(はず)。
いちいち書くのが面倒なので、次のようにマクロにしてしまうパターンが多い。

#define SIZE(a) sizeof(a)/sizeof(a[0]) // 配列名を引数として取る

ポインタを出さないと記事名詐欺になってしまうので使い道を考えてみる。

もちろん基本型とか構造体へのポインタを配列に入れて上のような使い方をするのも便利なのだが、あまり新しい感じはしない。
むしろ、このテクニックは関数ポインタと結びついたときに真価を発揮する。

いま、0以上の整数入力に対して、1個づつ異なった処理を行うコードを書いてみる。
またif文を使って普通に書いてみる。

void f0(int arg){ ... }
void f1(int arg){ ... }
void f2(int arg){ ... }
...
int main(){
int a = ?;//関数に渡す引数
int in; //入力を受け取る変数
...
if(in == 0){f0(a);}
else if(in == 1){f1(a);}
else if(in == 2){f2(a);}
...

同じような処理を何回も書くのは癪だ。

そこで、関数ポインタの出番だ。
同じ戻り値の型、引数の型・数を持つ関数なら(定義は個別にしなければならないが)これでまとめられる。

void f0(int arg){ ... }
void f1(int arg){ ... }
void f2(int arg){ ... }
...
int main(){
int a = ?;//関数に渡す引数
int in; //入力を受け取る変数
void (*f[])(int) = { f0, f1, f2, ... };
...
if( in >= 0 && in < sizeof(f)/sizeof(f[0]) ){
 f[in](a);
}
...

この記述の素敵なところは、何個も何個も同じような記述を繰り返さなくて済むことだ。
関数ポインタを使えば、引数を書く回数も減らせる。

書く回数を減らすことは、ミスを減らすことにつながる。

2016/2/3 コードに間違いがあったので修正

ポインタを学ぶ(3) ポインタの用途

C/C++ プログラミング

ポインタの実際の用途について考える。
とはいっても、主なのはこの3つくらいだろう。

関数から関数外の値を操作する

この用途でもっともよく知られているのは、swap関数の実装だ。
swap関数とは、2つの同じ型の引数を取り、その2つを交換する関数である。
(C++では<algorithm>(C++11からは<utility>)にテンプレート関数として用意されている)

今、int型引数のswap関数を実装したい。
ポインタを知らなければ次のように書こうとするだろう。

void swap(int a, int b){ // 間違った例
 int c = a;
 a = b;
 b = c;
}

この処理は、関数に固めないでmainにべた打ちする、とかならちゃんと動作する。
でも、関数にしてしまうとうまく動かない。
なぜなら、C/C++では(C++の参照を除いて)すべての引数は値渡しされるので、関数内で扱う引数は、与えられた引数のコピーに過ぎないからだ。

そこでポインタの登場だ。ポインタを使った正しい実装は次のようになる。

void swap(int* a, int* b){ // 正しい例
 int c = *a;
 *a = *b;
 *b = c;
}

引数でポインタを取るようになっている。これで引数の(コピーでない)実体にアクセスして変更することができる。
この関数を使うには、たとえばx,yを交換したいなら、swap(&x,&y)とすればいい。

ファイル操作をする時に、FILE型のポインタが必要なのも、これと同じ理由と推測できる。
構造体も値渡し(コピー)だから、操作にはポインタが必要なはずだ。

メモリ動的確保

ポインタがもっとも必要とされる場所と言えるだろう。

メモリ動的確保とは、言ってしまえば適当な要素数の配列を実行中に確保すること、と言えるだろう。
動的確保が普通の配列とどう違うかと言うと、要素数を実行中に決められることだ。

具体的なやり方は、ここには書かないことにしたい。

無駄な領域を使わないようにする

巨大な変数を引数としてとる関数を考える。
普通の値渡しで書くと、関数が実行される際には実引数がコピーされて、関数内で用いられる。
そうなると、同じデータが2つメモリ上にあることになり、無駄である。特に大きいものがコピーされる場合は。

そんなときに、ポインタを使う。
ポインタの大きさは、メモリアドレスの大きさ(32bitとか64bitがせいぜい)で、ほとんどの巨大な変数よりも小さい。

勉強不足がいよいよ誤魔化せないレベルになってきた。

2015/3/22 領域を節約する用途を追記(一応書いておいたほうがいい気がしたので)

ポインタを学ぶ(2) 特別なポインタ

C/C++ プログラミング

ポインタの実際の用途について学ぶ前に、少し特別なポインタについて学んでみる。

配列とポインタ

C/C++の配列はポインタによって実装されている。
たとえば、int a[20]を宣言したとする。
これが宣言されると、メモリ上にint型の領域が20個連続した領域に確保される。
連続した領域というのがミソで、これによって、先頭要素のアドレスさえわかればそのアドレスに加算していけば任意の要素にアクセスできる。
こういうわけで、配列先頭要素のアドレスは重要なので、配列名を単独で書くと配列先頭要素へのポインタになる。

この文法によって、配列要素へのアクセスにはいくつかの書き方がある。具体的にはこうだ。

// iは整数
*(a+i)
a[i]

正確には、上の書き方がむしろ本来の書き方で、下の書き方はその簡略版である、と言ったほうがいいかもしれない。

aとiは逆にしてもいいので、こんなのも同じ意味になる。

i[a]

競技プログラミングくらいでしか使わないだろうが。

関数ポインタ

実は、関数もメモリ上に場所をとるので、アドレスが存在する。
関数へのポインタは書き方がやや複雑だ。
次のような関数があるとする。

//Tは型名
T test(void)

別に引数がvoidである必要はない。説明を簡単にするためである。

この関数のアドレスを取得するには次のようにする。

test
// &test とも書ける

配列名と似たような感じで、関数名も関数へのポインタになるのだ。

アドレスがある以上、これを格納するポインタもある。
ある関数へのポインタの宣言は、次のように書ける。

戻り値の型(*ポインタ名)(引数の型);

さて、関数ポインタを宣言した。
これをどう使うかというと、結局のところは関数を呼び出すのに使うのだ。
関数ポインタに格納された関数を呼び出すには次のようにする。

(*ポインタ名)(引数);
// ポインタ名(引数); と書くこともできる


具体的に、関数ポインタtpを宣言して、test関数を格納して、呼び出してみる。

T(*tp)(void); //関数ポインタ tp の宣言
tp = test; // tp に test のアドレスを代入
(*tp)(); // tp() とも書ける。tpに格納されたアドレスにある関数を呼び出し

次回こそは、ポインタの実際の用途について学ぶ。

2014/3/14 関数ポインタから関数を呼び出す方法を追記
2016/2/3 解説を詳しくした

ポインタを学ぶ(1) ポインタの基本

C/C++ プログラミング

巷のC初心者に立ちはだかる最後にして最大の壁ともいわれる、ポインタと言う概念。
いったいポインタとは何なのか、何に使うのか、学んでいく。

基本

Cにおけるポインタとは

あるオブジェクト*1のメモリアドレスを値として持つ変数。
格納されているアドレスにあるオブジェクトを「指し示す」ことからポインタ(pointer)と呼ばれる。

ポインタを宣言するには、次のようにする。

型名 *ポインタ名

間接演算子*を使うことで、ポインタに格納されているアドレスにあるオブジェクトにアクセスできる。

*ポインタ名 //ポインタが指し示すオブジェクトにアクセス

メモリアドレスとは

あるオブジェクトがメモリ上のどこにあるかを指定するための一意な識別子。
オブジェクトのメモリアドレスを得たい場合は、アドレス演算子&を使う。

&オブジェクト名 //オブジェクトのメモリアドレスを得る

コードの例

#include <stdio.h>

int main(){
 int a = 10;        // int型変数 a を宣言し、初期値 10 で初期化
 int *p;            // int型に対するポインタ (int*型) の変数を宣言
 p = &a;            // p に a のアドレスを代入
 printf("%d\n",*p); // p の示すアドレスにある値を表示
}


このコードでは、(説明目的でなければ)ポインタを使う意味は微塵も無い。
ではポインタはどういう目的で使うのか。
それはまた次回。

2016/2/3 解説を文章中からコード内に移した

*1:変数とか関数のこと

最小二乗法(6) 線形近似における係数の推定誤差

物理・数学

前:最小二乗法(5) 線形近似 - equal_l2’s blog
次:最小二乗法(7) 項の欠けた線形近似 - equal_l2’s blog

\beta_0,\beta_1の推定誤差\sigma(\beta_0),\sigma(\beta_1)を求める。

「最小二乗法(2)」で示した通り、
\displaystyle \sigma(\beta_j)=\sigma_y \sqrt{\sum_{i=1}^{n}\left\{ \left( \frac{\partial \beta_j}{\partial y_i} \right)^2\right\}}
ただし、
\sigma_y=\displaystyle \sqrt{\frac{\sum_{i=1}^{n}\bigl(y_i- f(x_i)\bigr)^2}{n-d-1}} = \sqrt{\frac{\sum_{i=1}^{n}\bigl(y_i- f(x_i)\bigr)^2}{n-2}}

また、「最小二乗法(5)」で示した通り
\begin{eqnarray} \beta_0 &=& \frac{X_2 Y_0 - X_1 Y_1}{X_0 X_2 - {X_1}^2} \\ \beta_1 &=& \frac{X_0 Y_1 - X_1 Y_0}{X_0 X_2 - {X_1}^2} \end{eqnarray}

計算に入る前に、結論を言っておくと、
\begin{eqnarray} \sigma(\beta_0) &=& \sigma_y \sqrt{\frac{X_2}{X_0 X_2 - {X_1}^2}} \\ \sigma(\beta_1) &=& \sigma_y \sqrt{\frac{X_0}{X_0 X_2 - {X_1}^2}} \end{eqnarray}
である。

さて、計算に入ろう。

続きを読む

最小二乗法(5) 線形近似

物理・数学

前:最小二乗法(4) 多項式フィッティング - equal_l2’s blog
次:最小二乗法(6) 線形近似における係数の推定誤差 - equal_l2’s blog

\def \vb {\vec{\beta}} \def \fxb {f(x_i,\vb)}前回言った通り、d=1の場合について実際に解いてみる。

この場合のモデル関数は、あえて書くなら次のようになる。
\fxb= \beta_0 + \beta_1 x_i
もちろん、
\vb=\left(\begin{array}{c}{\beta_0 \\ \beta_1}\end{array}\right)
である。

このとき、行列形式の正規方程式は、
\begin{eqnarray} \left(\begin{array}{c}{Y_0 \\ Y_1}\end{array}\right) &=& \left( \begin{array}{cc} X_0& X_1 \\ X_1& X_2 \\ \end{array} \right) \left(\begin{array}{c}{\beta_0 \\ \beta_1}\end{array}\right) &=& \left( \begin{array}{cc} X_0& X_1 \\ X_1& X_2 \\ \end{array} \right) \vb \end{eqnarray}
ただし、
{\displaystyle X_m = \sum_{i=1}^{n}{x_i}^m,Y_m=\sum_{i=1}^{n}{x_i}^m y_i}
であった。

実際の係数がほしいだけなら、これに数値を入れて計算すればいい。

今回は、係数の誤差の式を求めるのに、各係数がX_i,Y_iにどのような依存をしているかを知る必要があるので、X_i,Y_iを使う形のまま解く。

次のような拡大係数行列
\left( \begin{array}{cc|c} X_0& X_1& Y_0 \\ X_1& X_2& Y_1 \\ \end{array} \right)
を行基本変形して\vbを求める。(掃き出し法)
(もちろん解さえ求まるなら方法は問わない。私は最初にやった時はクラメルの公式を使った)

結論を先に言えば、
\begin{eqnarray} \beta_0 &=& \frac{X_2 Y_0 - X_1 Y_1}{X_0 X_2 - {X_1}^2} \\ \beta_1 &=& \frac{X_0 Y_1 - X_1 Y_0}{X_0 X_2 - {X_1}^2} \end{eqnarray}
となる。

一応、以下に行基本変形の過程を逐一書いておく。わかる人には読むだけ時間の無駄だと思う。

続きを読む