2016-07-24
乱数とシミュレーション

作者: ryamada
発売日: 2016/07/22
メディア: Kindle版
この商品を含むブログ (1件) を見る
---
title: "RandomValues and Simulation"
author: "ryamada"
date: "2016年7月23日"
output: html_document
---

# Generate Random Values 乱数を作る

## Uniform distribution 一様分布

More numbers, more uniform. 数が増えると一様に近づく。
Generate many numbers and use a part of them for small sample size exapmles. 一度にたくさんの乱数を作っておいて、その一部を少サンプルの例として使う。

```{r}
n <- 10^6
r <- runif(n)
par(mfrow=c(2,3))
for(i in 1:6){
  tmp <- r[1:(10^i)]
  hist(tmp)
}
par(mfrow=c(1,1))
```

Sample mean converges to 0.5. 標本平均は0.5に近づく。

```{r}
n <- 10^4
r <- runif(n)
means <- rep(0,n)
for(i in 1:n){
  tmp <- r[1:i]
  means[i] <- mean(tmp)
}
plot(means,type="l",ylim=c(0,1))
```

Some functions work on vectors quicker than loop.
ベクトルに適用する関数には、ループを回すより速いものがある。

```{r}
# quicker
cmsm <- cumsum(r)
means.2 <- cmsm/(1:n)
plot(means.2,type="l",ylim=c(0,1))
```

## Normal distribution 正規分布
```{r}
n <- 10^4
hist(rnorm(n))
```

## Multidimensional 多次元

## (Square, Cube,...) vs. (Circle, Sphere,...) (正方形、立方体、…) 対 (円、球、…)

Square uniform distribution in 2D 2次元で正方領域一様分布。

Two-dimensional normal distribution. 2次元正規分布

```{r}
d <- 2
n <- 10^4
xy1 <- matrix(runif(n*d),ncol=d)
plot(xy1)

xy2 <- matrix(rnorm(n*d),ncol=d)
plot(xy2)
```

3D cubic. 3次元立方体。3D normal. 3次元正規分布
。
```{r setup}
library(knitr)
library(rgl)
knit_hooks$set(rgl = hook_rgl)
d <- 3
xy1 <- matrix(runif(n*d),ncol=d)
```

```{r, rgl=TRUE}
plot3d(xy1)
```
```{r, rgl=TRUE}
xy2 <- matrix(rnorm(n*d),ncol=d)
plot3d(xy2)
```

## Directionally Uniform 方向一様

Directional evaluation needs directinally uniform distribution.
方向についての評価では、方向に関する一様分布が必要になる。

Normal distribution is equal to all directions in any dimension.

正規分布は次元によらず、全方向について平等。

For 2D, uniform distribution of angle can be used.

2次元では、角度に関する一様分布も使える

```{r}
d <- 2
xy2 <- matrix(rnorm(n*d),ncol=d)
L <- sqrt(apply(xy2^2,1,sum))
xy2. <- xy2/L
plot(xy2.)
theta <- runif(n) * 2 * pi
xy3. <- cbind(cos(theta),sin(theta))
plot(xy3.)
```

3D. 3次元。

```{r, rgl=TRUE}
d <- 3
xy2 <- matrix(rnorm(n*d),ncol=d)
L <- sqrt(apply(xy2^2,1,sum))
xy2. <- xy2/L
plot3d(xy2.)
# ...
```

1D directionally uniform distribution...
1次元の方向一様分布とは？

```{r}
d <- 1
xy2 <- matrix(rnorm(n*d),ncol=d)
L <- sqrt(apply(xy2^2,1,sum))
xy2. <- xy2/L
plot(xy2.)
xy3. <- sample(c(-1,1),n,replace=TRUE)
plot(xy3.)
```

## Disk (circle inside) 円板

Disk is a peripheral circle with inside.
円板は円周と内部。

Uniform distribution of disk. 円板の一様分布。

Take random values only satisfying the rule. 条件に合う乱数のみを採用する。

```{r}
d <- 2
xy1 <- matrix(runif(n*d)*2-1,ncol=d)
L <- sqrt(apply(xy1^2,1,sum))
disk <- which(L<=1)
xy1. <- xy1[disk,]
plot(xy1.)
length(disk) / n
```

3D 3次元。

Fraction meeting the rule get lower.
採用割合が下がる。
```{r, rgl=TRUE}
d <- 3
xy1 <- matrix(runif(n*d)*2-1,ncol=d)
L <- sqrt(apply(xy1^2,1,sum))
disk <- which(L<=1)
xy1. <- xy1[disk,]
plot3d(xy1.)
length(disk) / n
```

10D's fraction is really low. 10次元の合格乱点割合は本当に低い。
```{r, rgl=TRUE}
d <- 10
xy1 <- matrix(runif(n*d)*2-1,ncol=d)
L <- sqrt(apply(xy1^2,1,sum))
disk <- which(L<=1)
xy1. <- xy1[disk,]
plot3d(xy1.)
length(disk)
length(disk) / n
```

## Discrete sampling 離散サンプリング

{0,1} even. {0,1}均等。
```{r}
n <- 10^4
s1 <- sample(0:1,n,replace=TRUE)
table(s1)
```

{0,1} with specific fraction.　割合を指定して{0,1}。

```{r}
n <- 10^4
p <- c(0.3,0.7)
s1 <- sample(0:1,n,replace=TRUE,prob=p)
table(s1)
```

n categories. カテゴリ数n。

```{r}
n <- 10^4
p <- c(0.5,0.2,0.1,0.1,0.1)
s1 <- sample(1:5,n,replace=TRUE,prob=p)
table(s1)
```

Number of samples per category can be generated with multinominal random generator. カテゴリ別のサンプル数のみを発生するなら多項分布乱数発生関数でできる。

```{r}
rmultinom(1,n,prob=p)
rmultinom(3,n,prob=p)
```

## Shuffle シャッフルする

sample() function can shuffle all, take out a part with or without replacement.

sample()関数は全要素の入れ替え、一部要素の取り出し(重複サンプリングについて、ありの場合となしの場合と)ができる。
```{r}
h <- c(1,4,5,7,9)
sample(h)
sample(h,replace=TRUE)
sample(h,3)
sample(h,10,replace=TRUE)
sample(h,4,replace=TRUE)
```

## Seeds　シード

Seeds control random value generators. シードによって乱数発生関数をコントロールする。

```{r}
runif(3)
runif(3)
set.seed(12345)
runif(3)
runif(3)
set.seed(12345)
runif(3)
runif(3)
set.seed(12345)
runif(3)
set.seed(12345)
runif(3)
runif(3)
set.seed(12345)
runif(6)

set.seed(12346)
runif(3)
```
## Random generating functions 乱数発生関数のいろいろ

```{r}
help(Distributions)
```
dxxxx: probability density function 確率密度
rxxxx: random value generator 乱数発生

# Simulation シミュレーション

## Damages pile up 傷害が積み重なる

Assume a cell will get damages periodically with probability p per a day.
When the number of damaging events reach q, the cell dies.
Generate distribution of life spans simulationally.

細胞はダメージを受けるが、そのダメージの生起確率は１日当たりpだと云う。ダメージイベントの回数がqになると細胞は死ぬという。
細胞の生存期間の分布をシミュレーションで作れ。

```{r}
n <- 10^4
p <- 0.05
q <- 20
days <- 10^3
x <- matrix(sample(0:1,n*days,replace=TRUE,prob=c(1-p,p)),ncol=n)
#image(x)
x.cumsum <- apply(x,2,cumsum)
matplot(x.cumsum[,1:100],type="l")
abline(h=q)
```
```{r}
plot(x.cumsum[,1])
abline(h=q)
which(x.cumsum[,1]==q)
min(which(x.cumsum[,1]==1))
res <- apply(x.cumsum,2,function(x){min(which(x==q))})
hist(res)
plot(sort(res))
```

## Differential equations with random errors ランダムエラー付きの微分方程式

$$
\frac{dx}{dt} = p x - q xy\\
\frac{dy}{dt} = r xy - s y
$$

```{r}
n.t <- 10000
p <-3
q <- 2
r <- 4
s <- 1

x.init <- c(1,1)
dt <- 10^(-3)

x <- matrix(0,n.t,2)
x[1,] <- x.init
for(i in 2:n.t){
  dxdt <- dt*(p *x[i-1,1] - q * x[i-1,1]*x[i-1,2])
  dydt <- dt*(r * x[i-1,1] * x[i-1,2] - s * x[i-1,2])
  x[i,1] <- x[i-1,1] + dxdt
  x[i,2] <- x[i-1,2] + dydt
}
matplot(x,type="l")
plot(x,type="l")
```

```{r}
x <- matrix(0,n.t,2)
x[1,] <- x.init
for(i in 2:n.t){
  dxdt <- dt*(p *x[i-1,1] - q * x[i-1,1]*x[i-1,2])
  dydt <- dt*(r * x[i-1,1] * x[i-1,2] - s * x[i-1,2])
  x[i,1] <- x[i-1,1] + dxdt*(1+rnorm(1)*1)
  x[i,2] <- x[i-1,2] + dydt*(1+rnorm(1)*1)
}
matplot(x,type="l")
plot(x,type="l")
```