---
title: "2次元ランダムウォーク"
author: '川久保圭吾'
date: "2015年10月12日"
output:
  html_document:
    number_sections: yes
    self_contained: no
    toc: yes
---

<link rel="stylesheet" href="https://rawgithub.com/brentertz/scianimator/master/assets/css/scianimator.css"/>
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.7.1/jquery.min.js"></script>
<script src="https://rawgithub.com/brentertz/scianimator/master/assets/js/jquery.scianimator.min.js"></script>

```{r setup, include=FALSE}
# アニメーション作成用のフック関数の指定
library(knitr)
opts_knit$set(animation.fun = hook_scianimator)
```

---

#ランダムウォークとは

確率変数
$$
X_1,X_2,\ldots,X_n \in \mathbf{R}^d \quad (n \in \mathbf{N})
$$
が独立同分布のとき，
$$
S_n = \sum_{i=1}^n X_i
$$
を*（d次）ランダムウォーク*という．

##例
$S_0 = O$ (原点) に置いた粒子が $i$ 回目の移動で $X_i \in \mathrm{R}^d$ と動いたときの位置が $S_n$であるとしたとき， $X_i$ が独立同分布ならば粒子の位置はランダムウォークである．

#ランダムウォークの再帰性

ランダムウォークについて，再帰的かどうか（初期位置に戻ってくるかどうか）という考えがある．
$n \to \infty$とするとき，初期位置に1回でも戻る確率が1のときに再帰的であると言う．
言い換えると，最初に初期位置に戻ってきた時刻を$T$とすると
$$P(T < \infty) = 1$$
である．

ここでは閾値$e$を設定し，ランダムウォークの最中に初期位置との距離が$e$未満になったら戻ったと判定する．
後で$\sqrt{\ }$を計算させる手間を省くため，あらかじめ2乗しておく．
```{r}
e <- (0.01)^2 # 戻ってきたと判定する距離
```

移動回数を$i$，試行回数を$j$とし
$$p_{i,j} := \frac{\hbox{i回のランダムウォークで戻ってきた回数}}{j}$$
を調べる．
今回のシミュレーションでは$i = imax = 10^4$，$1 \leq j \leq jmax = 10^5$とする．
```{r}
imax <- 10^4 # iの最大値
jmax <- 10^5 # jの最大値
```

```{r}
p <- rep(NA, jmax) # p_{imax,j}

check <- NULL # 戻ってきたかの判定（TRUE か FALSE）を格納する変数

MAX <- 10^6 # for 内の runif で一度に発生させる乱数の個数
once <- MAX %/% imax # 一度に行うランダムウォークの回数
times <- jmax %/% once # 一度に once 回行うランダムウォークを繰り返す回数

# S_1 から S_imax までの x 座標， y 座標，距離
x <- matrix(NA, imax, once)
y <- matrix(NA, imax, once)
s <- matrix(NA, imax, once)

for(i in 1:times){
    t <- matrix(2*pi*runif(MAX), imax, once) # i 行は i 回目に移動するときの角度

    # 列ごとに cos(t_i) と sin(t_i) の cumsum を求める
    x <- apply(cos(t), 2, cumsum)
    y <- apply(sin(t), 2, cumsum)

    s <- x^2+y^2 # 各距離を計算

    # apply(s < e, 2, any) は列ごとに，1回でも距離の2乗が e 以下になったかを論理値で返す
    # 1行 once 列で返ってくるので，それを1列に繋げてcheckに入れる
    check <- cbind(check, apply(s < e, 2, any))
}

p <- cumsum(check)/1:jmax # 各試行回数での計算結果を格納
```

ランダムウォークのある1回の結果
```{r random_walk-fig, fig.show="animate", interval=0.5, fig.width=4, fig.height=4}
for(n in seq(1, imax, imax %/% 30)){
    plot(x[1:n,1], y[1:n,1],
         xlim = c(min(x[,1]),max(x[,1])),
         ylim = c(min(y[,1]),max(y[,1])),
         type="l",
         xlab = "x",
         ylab = "y")
}
```


最大回数まで試行を行った結果は次のようになった．
```{r}
print(p[jmax])
```

順にプロットしていくと，最初は大きく振れるものの，徐々に1未満の値に収束していくように見られる．
```{r random_walk-fig2, fig.show="animate", interval=0.5, fig.width=4, fig.height=4}
for(n in seq(1, jmax, jmax %/% 30)){
    plot(p[1:n],
         xlim = c(1,jmax),
         ylim = c(0,max(p)),
         type="l",
         xlab = "Trial",
         ylab = "Probability")
}
```

