俺総研：So-netブログ

	ブログをはじめるログイン

2014年06月｜ 2014年07月｜2014年08月ブログトップ

二項分布から拡散方程式　 [Math]

二項分布では n ステップ後に位置 x にいる確率は p(x,n) で表しました。
この確率から差分方程式を作って拡散方程式が導かれます。せっかく p(x,n)
を前回計算したので、拡散方程式も見てみないとモッタイナイのでここに
書いておきます。なお、内容は検索で出てきた立命館大学の講義資料の
まんまです。

今回 x は変数変換後の変数として使うので、変換前は m と書いておきます。
なので、n ステップ後に位置 m にいる確率を p(m,n) と今回は表記します。

ステップ n+1 で m にいる確率 p(m,n+1)は以下のようになる。

$p(m,n+1)=\frac{1}{2}p(m+1,n)+\frac{1}{2}p(m-1,n)$
(ステップ n で m-1 にいる確率)*(右に移動する確率)と
(ステップ n で m+1 にいる確率)*(左に移動する確率の和。

これを変形して(こんな変形思いつかないけど、賢い人は変形できるんでしょうねぇ)
$p(m,n+1)-p(m,n)=\frac{1}{2}\left\{p(m+1,n)+p(m-1,n)-2p(m,n)\right\}$

変数変換で極限をとって連続にする。
$x=m\Delta x$
$t=n\Delta t$
$\Delta{x}\to{0},\Delta{t}\to{0}$

連続版の確率密度分布を f(x,t) とすると、区間 Δx の確率は次のように表せる。
$p(x,n)=f(x,t)\Delta{x}$

差分方程式を f(x,t) で表すと
$f(x,t+\Delta{t})=\frac{1}{2}\left\{ f(x+\Delta{x},t)+ f(x-\Delta{x},t)-2f(x,t)\right\}$
微分の形に変形して
$\frac{f(x,t+\Delta{t})-f(x,t)}{\Delta t}=\frac{{\Delta{x}}^2}{2\Delta{t}} \frac{ \frac{f(x+\Delta{x},t)-f(x,t)}{\Delta x} + \frac{f(x,t)-f(x-\Delta{x},t)}{\Delta{x}}}{\Delta x}$
こんな変形絶対できないんですけど。まぁ、天才アインシュタインが
考えてノーベル賞を取った理論だそうなので、できなくて当然。

$\frac{{\Delta{x}}^2}{2\Delta t}=D$ を保ったまま $\Delta x \to 0, \Delta t \ to 0$
そんなこと言われても、ああそうですかとしか。

最終的に得たのが次の式で、Dを拡散係数と呼ぶ。
$\frac{\partial f(x,t)}{\partial t}=D\frac{\partial^2 f(x,t)}{\partial x^2}$

元の講義ノートから変数の文字を変えたりして結局見にくくなった感が。。。
まぁこれを書いてわかった気になるのがひとつの目的なので。

2014-07-24 01:30 nice!(0) コメント(0) トラックバック(0)

二項分布でランダムウォーク　 [Math]

何かと出てくるランダウウォークですが、Scilab でシミュレーションして
理解した気になりたいところですが、こればかりは式の導出を避けては通れ
なさそうなので、導出をトレースしてわかった気になろうと思います。

原点からステップ毎に確率1/2で -1 or 1 を選んで移動する現象があるとします。
ステップ回数が多いとだんだん原点から離れていきます。この現象の何が大事
なのかというとこのサイトにとてもわかりやすくまとまっていて次の3つが
挙げられています。

1. n ステップ後に x にいる確率は正規分布になる
2. 分布の幅が√n で広がっていく
3. 原点にいる確率は 1/√n で減っていく

この3つを見ておきたいと思います。

式の導出は一応自分でも手計算で確認しつつ、こちらにまとめました。

nステップ後に位置mにいる確率
$p(m,n)=\frac{n!}{(\frac{n+m}{2})!(\frac{n-m}{2})!}(\frac{1}{2})^n$

式変形により正規分布の形に近似した式
$p(x,n)=\frac{1}{\sqrt{2\pi na^2}}exp(-\frac{x^2}{2na^2})$
a はステップ毎の移動距離。xは位置(x=a*m)。

平均0、分散μの正規分布
$p(x,n)=\frac{1}{\sqrt{2\pi \sigma^2}}exp(-\frac{x^2}{2\sigma^2})$

正規分布と比べると
$na^2=\sigma^2$
と対比できるので、「分布の幅が√n で広がっていく」のがわかります。

また、原点にいる確率は
$P(0,n)=\frac{1}{\sqrt{2\pi na^2}}$
これから「原点にいる確率は 1/√n で減っていく」というのも見て取れます。
なるほど。

2014-07-21 16:10 nice!(0) コメント(0) トラックバック(1)

スターリングの公式　 [Math]

最近個人的な理由で、ランダムウォークを理解するため、
二項分布から正規分布を求める式変形を理解しようとしています。
その中で出てくるスターリングの公式です。

ここで書くまでもなく、ググれば良質な情報がいくらでも出てきますが、
理解した気になるために、ここに書いておきます。

$n! \approx n^n e^{-x} \sqrt{2\pi n}$

両辺の対数をとる

$\ln{n!} \approx n\ln{n} -n + \frac{1}{2} \ln{{2\pi} + \frac{1}{2} \ln{n}$

2014-07-20 00:38 nice!(0) コメント(0) トラックバック(0)

log(1+x)をマクローリン展開　 [Math]

二項分布から正規分布を求めようとしています。
その中で log(1+x) のマクローリン展開が出てくるのでメモっておきます。

マクローリン展開はテイラー展開の x=0 の周りで展開するバージョン。

テイラー展開
$\Sigma_{i=0}^{\infty}\frac{f^{(n)}(a)}{n!} (x-a)^n$

マクローリン展開は a=0 なので
$\Sigma_{i=0}^{\infty}\frac{f^{(n)}(0)}{n!} x^n$

n=3 以降は O(3) とかやって捨てるので n=2 までが必要です。
そのためには x=0 のときの log(1+x) の微分と二回微分が出てくるので計算しておきます。

log(1+x)の微分
$\frac{d}{dx}{\log(1+x)}=\frac{1}{1+x}$
x=0 のときは 1

log(1+x)の二回微分
$\frac{d^2}{dx^2}{\log(1+x)}=-\frac{1}{(1+x)^2}$
x=0 のときは -1

あと log(1+x) の x=0 のときのゼロ階微分は 0。
これらを使って、log(1+x) をマクローリン展開

$\log(1+x) = \Sigma_{i=0}^{\infty}\frac{\log(1+x)^{(n)}(0)}{n!} x^n\\= 0 + \frac{1}{1!}x + \frac{-1}{2!}x^2 + O(3) \\ \approx x-\frac{x^2}{2}$

3乗の項以降は O(3) と書いて、最終的には捨てました。

2014-07-20 00:24 nice!(0) コメント(0) トラックバック(0)

Scilab でコイン投げの二項分布　 [Math]

Scialb でシミュレートしてわかった気になるシリーズ第三弾。
コイン投げを、n 回投げて、 x 回表が出る確率の分布が二項分布。

1.n回コインを投げる
2.表が出た回数を数える
....というのを1000回繰り返して分布を描く。

n=5

n=100

n が大きくなるにつれて、平均=n*p、分散=n*p(1-p) ただしp=1/2 の正規分布に近づいている。
なんか嘘っぽいけど。

Scilab のコード。Github にもおいてあります。

clear;

// sampling num
samplenum=100
// number of repetition
repeat=1000;

// coin toss
R=grand(repeat,samplenum,'uin',0,1);

scf(1);
sums=sum(R,'c');
histplot([0:1:samplenum],sums)

n=samplenum;
p=1/2;

avg=n*p;
s=sqrt(n*p*(1-p))
x=[0:0.1:samplenum];
plot(x,(1/(s*sqrt(2*%pi)))*exp(-((avg-x).^2)./(2*s^2)));

2014-07-18 02:05 nice!(0) コメント(0) トラックバック(0)

コイン投げで中心極限定理　 [Math]

scilab でシミュレーションしてわかった気になるシリーズ第二弾です。

コイン投げで表=1、裏=0とすると、平均値は0.5、分散は0.25。コイン投げというとなぜか n 回投げて x 回表が出る確率が二項分布になる話が多いですが、今回は平均0.5で分散が0.25/√n になるのをサクッと scilab で見てみることにします。

1.n回コインを投げる
2.平均値をプロット
....というのを100回繰り返す。

n を大きくしていくとばらつき(分散)が小さくなって、収束していく

n=10

n=100

n=500

ほんとだ、正規分布になっていて、 n が増えるほど分散が小さくなっている。

Scilab のコード。Github にもおいてあります。

clear;

// sampling num
samplenum=500;
// number of repetition
repeat=100;

// coin toss
R=grand(samplenum,repeat,'uin',0,1);
Avgs=mean(R,'r');
histplot([0:0.01:1],Avgs)

avg=1/2;
s=sqrt((1/4)/samplenum);
x=[0:0.01:1];
plot(x,(1/(s*sqrt(2*%pi)))*exp(-((avg-x).^2)./(2*s^2)));

2014-07-18 01:10 nice!(0) コメント(0) トラックバック(0)

Scilabで中心極限定理　 [Math]

中心極限定理という統計学の定理がありまして、何度聞いても
わかったような、わからないような。なので Scilab でサクっと
シミュレーションして、わかった気になろうという思いです。

標本平均は一定数の観測値を集めて平均したもの。サイコロなら100回
振って出た数の平均値のこと。

この平均値をいくつも集めてると正規分布になるという。
もとの現象の分布は関係ない。サイコロはどれが出るかは1/6で
一定に分布している。でも標本平均を何回もとると正規分布になる。

真の平均値μ、分散σ^2という現象を観測していたなら、
n個の標本平均をとると、標本平均の分布は、平均μ、分散はσ^2/√n
になるという。

サイコロの場合は平均=3.5、分散=2.92。これを使って Scilab で見てみます。

1000回サイコロを振ったヒストグラムをプロット。

それを10回繰り返し、1回あたりの平均を取った(10このサイコロを1000回振って、
各サイコロの平均を取ったのと同じ)。

繰り返しを200回にしてみた。
ついでに平均=3.5、分散=√2.92/200の正規分布を描いた。

あー、なってますね。正規分布に。

グラフ作成に使った scilab のコード

clc;

// number of repetition
repeat=10;
// draw dice 10 times. and repeat it.
X=grand(1000,repeat,'uin',1,6);
// avarage of repetition
R=mean(X,'c');
// plot histgram
histplot([0:0.1:6],R);
title="Frequency distribution of dice face with " +string(repeat)+" roll(s)";
xtitle(title,"dice face","Frequency");

avg=3.5
s=sqrt(2.92/repeat);
x=[0:0.01:6];
plot(x,(1/(s*sqrt(2*%pi)))*exp(-((avg-x).^2)/(2*s^2)));