2014年10月12日日曜日

数値計算の常識

少し前に読んだ本。

興味のある章だけ飛ばし読みしたんですが、丸め誤差と桁落ち、刻み幅と数列の収束の話は読んでよかったと思う。どう考えても常識として知っておくべき。

桁落ちについてはこれまで意識してなかったので、本の中にもあった刻み幅と差分の関係を調べてみると確かにその通りでした。
y=sin(x)のx=1での差分と真値cos(1)の差を出力。 (これだとcos(1)の丸め誤差も気になる?)
刻み幅は、最も小さい点ではマシンイプシロン。そこからは順に10倍しています。

やはり刻み幅は小さければいいというものではなく、ある程度小さくなると桁落ちの影響が相対的に大きくなっていくのが見て取れます。

前進差分と中心差分で二通りやってますが、中心差分の方が精度良いというのも教科書通り。

 しかし、最初計算してみたときはx=0で差分を計算してみたのですが、このときは差分がマシンイプシロンのときも正しい値を返していました。
どうやら桁落ちの影響は常にあるわけではないようで、テストの値をきちんと用意してやらないと桁落ちの影響を見逃してしまいそうだと思いましたまる

2014年9月27日土曜日

テンプレート引数を式パラメータにしたとき

TMPの練習にと、FizzBuzzしたり素数並べたりしてるとこで

template<int N, int D = N - 1> struct Test

っていう形に出会った。パッと見第2引数が常にN - 1になるかと思ったが、実はそうではなかった。

テンプレートの第2引数に何も指定されていなければD=N - 1になり、何か引数が与えらていればDはその値になるようだ。

例えば以下のコードの出力は
9
5
10
 になる。

#include <iostream>

template<int N, int D = N - 1> struct Test{
    Test(){
        std::cout << D << std::endl;
    }
};

template<int N> struct Test<N, 1>{
    Test(){
        std::cout << N << std::endl;
    }
};

int main()
{
    Test<10>();
    Test<10, 5>();
    Test<10, 1>();
}


ちなみに3つ目は部分特殊化

2014年9月8日月曜日

テンプレートメタプログラミングが難しい

Effective C++を読んでいたら最後の法でテンプレートメタプログラミング(TMP)というのが出てきました。

以下の階乗を計算するプログラムがTMPの初歩的なプログラムらしいです。

#include <iostream>

template<unsigned n>
struct  Factorial {
    enum { value = n * Factorial<n - 1>::value };
};

template<>
struct  Factorial<0> {
    enum { value = 1 };
};

int main()
{
    std::cout << Factorial<5>::value << std::endl;
    std::cout << Factorial<10>::value << std::endl;
    getchar();

}

このコードを見てよくわからなかったのは

・enumってなんか意味あるの?
・Factorialっていつ計算されるの?

でした。

enumってなんか意味あるの? 

enumハックというやつらしいです。コンパイル時にクラス内で定数が必要になるときに使うそうです。

static const unsigned value = n * Factorial<n - 1>::value;

と書いても動きましたが、こっちはコンパイラによってはエラーになるらしい。

Factorialっていつ計算されるの?

コンパイル時にされる。普通に本文に書いてあったのですが、文章ってあまり頭に入ってこないですね。コード見たときになんやこれ、って思って悩んでしまいました。
TMPでは、このような静的な計算がコンパイル時に行われるため、計算が高速化するらしいです。


ほかにTMPの効用としては

エラーがコンパイル時にみつかる

たとえば Factorial<20> とかにすると値が大きくなりすぎてエラーをはく

TMPはチューリング完全であることが証明されている

要はどんな計算でもTMPで書けることが保障されているらしい・・・



しかし、実際どうやってTMPで書いていったらいいかわからんです


2014年8月3日日曜日

プログラミング

4月から会社でプログラミングの研修を受けてました。

そんで一応動くコードが書けるようになってきたところ、はてぶの

クイズに答えるだけで「プログラミングの勉強」ができるWebサービスたちにハマる! 

 という記事を見つけました。

自分がどのくらいコードを書けるようになったのか知るにはちょうどいいかなーと思い、CodeEvalとpaizaをやってみました。

CodeEvalの方は、英語なのが若干めんどいですが、問題数がたくさんあってこれやってたらいつまでも時間を潰せそう。

paizaの方は問題を解くのに時間制限があって、しかも一回しか挑戦できないからなんかわからんけど緊張感がありました。なにやら書いたコードを企業に送って転職活動をするというのが主目的のサイトらしい。どうりで問題をブログとかに書くなよって警告されるわけです。

実際に問題やってみてわかったことは、プログラミングのスキルというより、入出力のスキルが不足しまくってるということでした。

問題を解く部分はわりとすぐできるのに、入力を受け取ることに時間かかりまくってました。

C++で書いてたのですが、スペース区切りの文字列を受け取るのができなくて発狂しそうでした。

ググってみたところgetlilne関数使えばいいよー的な記事があったので、getline(cin, str)って書いたら、そんな関数ねーよと言われてしまいました。ちゃんとヘッダもインクルードしたはずなのに。なんやねん

そんなこんなで入出力の部分を作ることができたら、そのあとはすんなり作ることができました。算数とか得意だったからそっちはいけるらしい。

 難しい問題も少し見てみたところ、問題になりそうなのはやはり入出力でした。標準入出力ができるようになったら難しい問題に挑戦してみたいと思います。

 

2014年5月20日火曜日

制御変量法のそれらしい解釈

制御変量法を用いることでなぜ分散が減少するのかについてそれっぽい解釈を思いついた

たとえばモンテカルロシミュレーションで求めた$Y$には真の値からの誤差$\delta _Y$が存在する。
また$\hat{\theta _c}=Y+c(Z-E(Z)) $の$Z$も真の値$E(Z)$からのずれ(誤差)$\delta _Z$が存在する。
これらを前記事(1)に代入すると

\begin{equation}
\hat{\theta _c}=Y+c(Z-E(Z))\\
=Y_{true}+\delta _Y +c(Z_{true}+\delta_Z-E(Z))
\end{equation}

$Z_{true}=E(Z)$とすると

\begin{equation}
\hat{\theta _c}=Y_{true}+\delta _Y +c\delta_Z
\end{equation}

となる。$\delta_Y$と$\delta_Z$はシミュレーションによって生じた真の値からの誤差であるが、YとZの間になんらかの相関があればこれら$\delta_Y$、$\delta_Z$の間にもそれに近い相関があると考えられる。

cの値を適切に選ぶことで$\delta _Y +c\delta_Z$の分散を小さくしようというのが制御変量法の発想だと思う。ふたつの相関する乱数である$\delta _Y$と$\delta_Z$に対して、cを前記事のように$c^*=\frac{\rm{cov}(Y,Z)}{\rm{var}(Z)}$と選ぶと分散を圧搾できると考えられる。


これと同じような考えで、モンテカルロシミュレーションである量$f(x)$の1階微分を求めるとき

\begin{equation}
\frac{df}{dx}=\rm{lim}_{\Delta x\rightarrow \infty}\frac{f(x+\Delta x)-f(x)}{\Delta x}
\end{equation}
の式を用いるが、$f(x+\Delta x)$と$f(x)$を求めるとき、同じ乱数列を用いるほうがよいのか、異なる乱数列を用いるほうがよいのかというお話がある。

さきほどの考えで行くと、$f(x+\Delta x)$と$f(x)$のどちらも真の値から誤差があるが、どちらの計算にも同じ乱数列を用いると、それらの誤差は$\Delta x$による差しか生じないためかなり相関していると考えられる。しかし異なる乱数列を用いるとこの相関はなくなってしまう。

 $\frac{df}{dx}$を求めるときには$f(x+\Delta x)$と$f(x)$の差をとるため、これらの誤差に相関があれば、差をとったときに誤差も相殺してしまうことができる。このため、$\frac{df}{dx}$の値の分散は同じ乱数列を用いたときのほうが異なる乱数列を用いたときより小さくなるだろう。

2014年3月29日土曜日

制御変量法

制御変量法について日本語で検索してみてもwikiで
"モンテカルロ法における分散減少法の1つ. 推定したい特性値の不偏推定量に対して, これと相関があって期待値がわかっている確率変数のことを制御変量という. 不偏推定量と制御変量の1次結合をうまく作れば,推定の分散を減らすことができる. "
程度しかわかりませんでした。

英語で検索してみると解説してくれるPDFがすぐに見つかって、読んでみると思ったより簡単そうだったので適当にメモ。

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

$\theta :=E(Y)$を決定したい。ここでYはある試行の結果であるとする。
また、Zを同じ試行から得られる量で、その期待値E(Z)は既知であるとする。
Zのことを制御変量(control variate)と呼ぶ。

これらについての2通りの不偏推定量
\begin{equation}
\hat{\theta }=Y \label{simple}
\end{equation}
\begin{equation}
\hat{\theta _c}=Y+c(Z-E(Z)) \label{contvar}
\end{equation}
を考える。ここで明らかに$E(\hat{\theta _c})=\theta$である。

考えるべきは$\hat{\theta}$の分散と$\hat{\theta_c}$の分散の大小関係である。分散が小さいほど、$\theta$についてより精度のよい推定値であるということになる。
\begin{equation}
{\rm var}(\hat{\theta_c})={\rm var}(Y)+2c{\rm cov}(Y,Z)+c^2{\rm var}(Z) \label{thetac}
\end{equation}
\begin{equation}
{\rm var} (\hat \theta )={\rm var}(Y)
\end{equation}
式\ref{thetac}をcについての二次式と考えると、$c=-\frac{{\rm cov}(Y,Z)}{{\rm var}(Z)}\equiv c^*$のとき${\rm var}(\hat{\theta_c})$は最小値
\begin{equation}
{\rm var}(\hat{\theta_{c^*}})={\rm var}(\hat \theta )-\frac{{\rm cov}(Y,Z)^2}{{\rm var}(Z)}
\end{equation}
をとることが分かる。これは${\rm cov}(Y,Z)\not =0$である限り${\rm var}(\hat{\theta})$より小さい。

したがって、$\theta :=E(Y)$を推定するときには、推定値として式\ref{simple}を用いるより式\ref{contvar}を用いたほうが得られる推定値の分散が小さくなる。分散の減少量は${\rm cov}(Y,Z)$の二乗に比例するが、共分散${\rm cov}(Y,Z)$はYとZの相関の大きさを表すので、相関が大きいほど、つまり制御変量Zによって得られるYについての情報量が多いほど精度よく$\theta $を推定できる。と感覚的にわかる。しかし${\rm cov}(Y,Z)$と${\rm var}(Z)$は独立でないために、${\rm var}(\hat{\theta_c})$を最小にするZは${\rm cov}(Y,Z)=\pm 1$とは限らないようだが、今はまあいいかと思ったので詳しくは調べていない。

実装に当たっては、$c^*$を使うために${\rm cov}(Y,Z)$と${\rm var}(Z)$を知る必要がある。
これらは既知とは限らないが、例えばn回目の試行なら、n-1回目までの試行で得られた値とE(Z)から推定した値を用いることで代用するらしい。

2014年3月17日月曜日

上極限と下極限

上極限と下極限について定義の式を見ただけでは全然意味がわからなかったです。いろいろ調べてみたんですが、結局ルベーグ積分30講の解説を読んで

xが$A_n$の下極限に含まれる$\Leftrightarrow $ある$n_0$が存在して$n \geq n_0 \Rightarrow x\in A_n$

xが$A_n$の上極限に含まれる$\Leftrightarrow$どんな$N>0$に対しても$n>N$かつ$x\in A_n$を満たすnが存在する

ということだと理解しました。下極限のほうがわかりやすかったので下極限、上極限の順になってます。

定義式を最初に見たとき変な形の式だなと思ったので、式的な意味でも理解しておくべきかと思い式を眺めていたら上極限と下極限を導入した経緯について自分なりに納得できたので以下にメモ。


以下は私が勝手に考えただけなので間違ってるかもしれません
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

集合列$\{A_n\}$、$A_1, A_2, A_3, ...A_n, ...$

について、この集合列の極限$\lim_{n \to \infty}A_n$を考える。

$A_n$は任意の集合列であるので、そもそも極限が存在することも定かではない。

そこで、先ずは極限が存在することがわかっている単調増加/単調減少の数列を出発点とする。

$B_n$を単調増加の数列、$C_n$を単調減少の数列とすると、これらの極限は
\begin{eqnarray}
\lim_{n \to \infty}B_n=\bigcup_{n=1}^{\infty}B_n
\\
\lim_{n \to \infty}C_n=\bigcap_{n=1}^{\infty}C_n
\end{eqnarray}
となる。このように単調増加/単調減少の数列ならば極限を考えることができるので、$\{A_n\}$をもとに単調増加/単調減少の数列を作り、それらの極限を考えることにする。

$\{A_n\}$をもとに単調増加の数列を作る。
\begin{equation}
B_k=\bigcap_{n=k}^{\infty}A_n
\end{equation}
は単調増加の数列である。$B_k$は$\{A_n\}$をn=kから無限大の積集合なので、kが大きくなるほど積をとる対象が少なくなるため、$k\leq l\Rightarrow B_k \subset B_l$となるのである。

$\{B_k\}$の極限を考える。
$B_k$は単調増加の集合列なので上の式より
\begin{equation}
\lim_{k \to \infty}B_k=\bigcup_{k=1}^{\infty}B_k
\end{equation}
である。この式は
\begin{equation}
\lim_{k \to \infty}B_k=\bigcup_{k=1}^{\infty}B_k=\bigcup_{k=1}^{\infty}\bigcap_{n=k}^{\infty}A_n
\end{equation}
となる。この、$A_n$をもとして作った単調増加の集合列$B_k$の極限を$A_n$の下極限と定義する。
\begin{equation}
\lim_{k \to \infty}B_k=\bigcup_{k=1}^{\infty}B_k=\bigcup_{k=1}^{\infty}\bigcap_{n=k}^{\infty}A_n\equiv  \varliminf_{n \to \infty} A_n
\end{equation}
下極限は以上のような手順で作ったため、$A_n$がどのような集合列であっても下極限は存在する。

単調増加列の場合と同様に単調減少列
\begin{equation}
C_k=\bigcup_{n=k}^{\infty}A_n
\end{equation}
を作り、その極限を上極限と定義する。
\begin{equation}
\lim_{k \to \infty}C_k=\bigcap_{k=1}^{\infty}C_k=\bigcap_{k=1}^{\infty}\bigcup_{n=k}^{\infty}A_n\equiv  \varlimsup_{n \to \infty} A_n
\end{equation}
こちらも$A_n$がどのような集合列であっても存在する。

これら$A_n$をもとに作った2種類の極限が一致するとき、それを$A_n$の極限と呼ぶ。
\begin{equation}
 \varlimsup_{n \to \infty} A_n=\varliminf_{n \to \infty} A_n \equiv \lim_{n \to \infty}A_n
\end{equation}

ルベーグ積分30講 (数学30講シリーズ)


読みました。

アマゾンで、これは読み物だから電車の中でもすらすら読めちゃうぜ~、というレビューを信じて購入。最初は読みやすかったけれど、後半はわかったような気はするけど実は全然わかっていないという学部時代の懐かしい感覚を思い出し、やはり数学の本だなーと感じました。
あとルベーグ積分の応用については物足りないどころか、読み終わったとき「で、何がしたかったんだっけ?」となってしまいました。

真面目に理解しようとしたら少なくともあと2回は読まないといけない程度の理解度ですが、ルベーグ積分の雰囲気だけは感じられたのでひとまず読了ということにします。あと2回読むには退屈しそう

2014年3月10日月曜日

ブログ作成

理系大学院修士を卒業した人が読んだ本とかについて書く予定です。

明日から旅行へ行くので帰ってきたら更新してみようと思います。