VBAユーザがPython・Rを使ってみた：乱数

# はじめに

機械学習の勉強を始めたVBAユーザです。
備忘録としてPython・Rの文法をVBAと比較しながらまとめています。

今回は、乱数についてまとめます。
確率分布については、[別記事](https://qiita.com/swathci/private/7d528fcb34a2e31a8127)にしました。

**目次**

- [乱数](#乱数)
   - [連続確率分布に従う乱数](#連続確率分布に従う乱数)
      - [連続確率分布に従う乱数の発生](#連続確率分布に従う乱数の発生)
      - [連続確率分布の例](#連続確率分布の例)
   - [離散確率分布に従う乱数](#離散確率分布に従う乱数)
      - [離散確率分布に従う乱数の発生](#離散確率分布に従う乱数の発生)
      - [離散確率分布の例](#離散確率分布の例)
- [まとめ](#まとめ)
  - [一覧](#一覧)
  - [プログラム全体](#プログラム全体)


# 乱数
今回は、乱数についてまとめます。

## 連続確率分布に従う乱数
まず、[連続確率分布](https://ja.wikipedia.org/wiki/連続確率分布)（continuous probability distribution）に従う乱数の発生から説明します。

### 連続確率分布に従う乱数の発生

例として、正規分布に従う乱数で説明します。

**Python**
Pythonでは、科学技術計算ライブラリの [SciPy](https://docs.scipy.org/doc/scipy/reference/) の統計関係のサブパッケージ [scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html) を使います。scipy.stats から使用する関数をインポートしておきます。（なお、標準ライブラリの[random](https://docs.python.org/ja/3/library/random.html)モジュールや、Numpyの[numpy.random](https://docs.scipy.org/doc/numpy-1.12.0/reference/routines.random.html)モジュールを使っても乱数を発生できます。これについては下参照。）
計算に使うので、Numpy もインポートしておきます。また、ヒストグラムを描画するために [Matplotlib](https://matplotlib.org/index.html) の [matplotlib.pyplot](https://matplotlib.org/api/_as_gen/matplotlib.pyplot.html?highlight=pyplot#module-matplotlib.pyplot) もインポートしておきます。

```Python:Python3
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

m = 0
s = 1
r = norm.rvs(m, s, size=1000, random_state=0)
plt.figure()
plt.hist(r, bins=50)
plt.show()
```

![histgram_norm1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/7898ae5b-257a-651a-aa71-7ec927ffb503.png)

乱数の発生はこれでできていますが、ヒストグラムの描き方について補足します。
`plt.hist(r)`でベクトル `r` のヒストグラムが描画できます。`bins=50` でビンの数をしています。
これに、正規分布の[確率密度関数](https://qiita.com/swathci/private/7d528fcb34a2e31a8127#確率密度関数)（pdf）を重ねて描画してみます。確率密度関数は全体が1になるように正規化されていますので、ヒストグラムも全体が1になるように正規化しておきます。`plt.hist`の引数で`density=True`と指定することで、全体を1に正規化したヒストグラムにできます。

```Python:Python3
m = 0
s = 1
r = norm.rvs(m, s, size=1000, random_state=0)
plt.figure()
plt.hist(r, bins=50, density=True, label='histgram') # 全体を1に正規化
x = np.arange(-5, 5, 0.1)
y = norm.pdf(x, m, s)
plt.plot(x, y, label='pdf')
plt.legend()
plt.show()
```

![histgram_norm2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/be14662b-ba4a-6bf9-dca0-8f077b3de9c7.png)


なお、Pythonでは他にも、Python標準ライブラリの[random](https://docs.python.org/ja/3/library/random.html)モジュールや、Numpyの[numpy.random](https://docs.scipy.org/doc/numpy-1.12.0/reference/routines.random.html)モジュールを使っても乱数を発生できます。
標準ライブラリの`random`モジュールでは1個の乱数を発生します。
Numpyの`numpy.random`モジュールでは、数ベクトルの形だけでなく配列の形の乱数も発生できます。
詳しくは下の[一覧表](#一覧)をご覧ください。

```Python:Python3_random
# random
import random
random.seed(0) # 乱数のシードの固定
random.normalvariate(0,1)

# 正規分布のヒストグラム
m = 50
s = 10
random.seed(0)
r = list()
for i in range(1000):
    r.append(random.normalvariate(m, s))
plt.hist(r, bins=50)
plt.show()

random.seed(0)
r = [random.normalvariate(m, s) for i in range(1000)] # リスト内包表記
plt.hist(r, bins=50)
plt.show()

r = np.zeros(1000)
random.seed(0)
for i in range(1000):
    r[i] = random.normalvariate(m, s)
plt.hist(r, bins=50)
plt.show()
```

```Python:Python3_numpy.random
# numpy.random
import numpy as np
np.random.seed(0) # 乱数のシードの固定
np.random.normal(m, s, 10)      # サイズ 10
np.random.normal(m, s, (10,10)) # サイズ (10,10)
np.random.normal(m, s, (2,3,4)) # サイズ (2,3,4)

# 正規分布のヒストグラム
m = 50
s = 10
np.random.seed(0)
r = np.random.normal(m, s, 1000)
plt.hist(r, bins=50)
plt.show()
```


**R**
Rでは、特に何もしなくても標準パッケージで確率分布が扱えます。
ヒストグラムは`hist()`で描画できます。

```R:R
m <- 0
s <- 1
set.seed(0) # 乱数のシードをセット
r <- rnorm(n=1000, mean=m, sd=s)
hist(r)
hist(r, breaks=50)
```
これに確率密度関数を重ねて描画するには、`par(new=T)`とした後に`plot()`します。Rでは関数の描画は`curve()`でもできます。

```R:R
m <- 0
s <- 1
set.seed(0)
r <- rnorm(n=1000, mean=m, sd=s)
xrange <- c(-4,4)
yrange <- c(0,0.6)
hist(r, breaks=50, prob=T, xlim=xrange, ylim=yrange, ann=F)
x <- seq(-5, 5, by=0.1)
y <- dnorm(x, mean=0, sd=1)
par(new=T)
plot(x, y, xlim=xrange, ylim=yrange, xlab="", ylab="", col="red", type="l")
par(new=T)
curve(dnorm, from=-4, to=4, ylim=yrange, col="blue", ylab="")
```
注意）Rでグラフを描画するには、標準パッケージでもできますが、[`ggplot2`](https://ggplot2.tidyverse.org/)を使った方が「一貫性のある文法で合理的に」描けます。（[こちら](https://heavywatal.github.io/rstats/ggplot2.html)の記事参照。）


**VBA**
VBAでは、0以上1未満の一様乱数なら[Rnd 関数](https://docs.microsoft.com/ja-jp/office/vba/language/reference/user-interface-help/rnd-function)で発生できます。引数の値によって乱数の生成方法が変わります（[こちら](https://docs.microsoft.com/ja-jp/office/vba/language/reference/user-interface-help/rnd-function#return-values)）。
[Randomize ステートメント](https://docs.microsoft.com/ja-jp/office/vba/language/reference/user-interface-help/randomize-statement)で、Rnd 関数の乱数ジェネレーターを初期化します。引数に数値を指定するとその数値を新しいシード値として使用され、引数を省略するとシステムタイマーから返される値が新しいシード値として使用されます。
正規乱数など、他の分布に従う乱数については、下の例を参照。

```VBA:VBA
Randomize 1
For i = 1 To 1000
    Debug.Print Rnd
Next i
```

数値の引数を指定した Randomize を使用する直前に、負の引数を指定した Rnd を呼び出しておくことで、同じ乱数のシーケンスを発生させることができます。

```VBA:VBA
Dim n As Integer
Dim r() As Single
n = 1000
ReDim r(1 To n)
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = Rnd
    Debug.Print i, r(i)
Next i
' 1             0.3335753
' 2             6.816387E-02
' 3             0.5938293
' 4             0.7660395
' 5             0.1892894
' 6             0.5373986
' 7             0.3269944
' 8             0.393937
' 9             7.341915E-02
' 10            0.8315425
```

### 連続確率分布の例

**Python**

```Python:Python3
from scipy.stats import uniform, norm, chi2, t, f, beta, gamma, lognorm, expon
```

- [一様分布](https://ja.wikipedia.org/wiki/一様分布) $U(a,b)$

```Python:Python3
# 一様分布 U(a,b)
a = 10
b = 20
plt.figure()
r = uniform.rvs(loc=a, scale=b-a, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(a, b, 0.1)
y = uniform.pdf(x, loc=a, scale=b-a)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = uniform.stats(loc=a, scale=b-a, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_uniform.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/40d8d626-d6ed-3b9c-8b38-c6931f10025d.png)


- [正規分布](https://ja.wikipedia.org/wiki/正規分布) $N(\mu,\sigma^2)$

```Python:Python3
# 正規分布 N(m,s)
m = 0
s = 1
plt.figure()
r = norm.rvs(loc=m, scale=s, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(-5, 5, 0.1)
y = norm.pdf(x, loc=m, scale=s)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = norm.stats(loc=m, scale=s, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_norm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/fbcbe168-a080-1293-c5c7-7f3444e68130.png)


- [t分布](https://ja.wikipedia.org/wiki/T分布) $t(n)$

```Python:Python3
# t分布 t(n)
n = 10
plt.figure()
r = t.rvs(df=n, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(-5, 5, 0.1)
y = t.pdf(x, df=n)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = t.stats(df=n, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_t.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/ce85cc31-40c7-0bd9-be63-52bc506ea3ab.png)


- [カイ二乗分布](https://ja.wikipedia.org/wiki/カイ二乗分布) $\chi^2(n)$

```Python:Python3
# カイ二乗分布 X2(n)
n = 10
plt.figure()
r = chi2.rvs(df=n, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 30, 0.1)
y = chi2.pdf(x, df=n)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = chi2.stats(df=n, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_chi2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/867f9a80-c86c-22cc-e87a-f610e25075c8.png)


- [F分布](https://ja.wikipedia.org/wiki/F分布) $F(n,m)$

```Python:Python3
# F分布 F(n,m)
n = 5
m = 20
plt.figure()
r = f.rvs(dfn=n, dfd=m, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 5, 0.1)
y = f.pdf(x, dfn=n, dfd=m)   # numerator(分子), denominator(分母)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = f.stats(dfn=n, dfd=m, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_f.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/7ec93040-479b-dc2f-c082-dbbf3d64f11e.png)


- [ベータ分布](https://ja.wikipedia.org/wiki/ベータ分布) $Be(a,b)$

```Python:Python3
# ベータ分布 Be(a,b)
a = 2
b = 5
plt.figure()
r = beta.rvs(a=a, b=b, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 1, 0.01)
y = beta.pdf(x, a=a, b=b)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = beta.stats(a=a, b=b, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_beta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/5e3c4a03-ae12-b0b0-730d-45845cdd734a.png)


- [ガンマ分布](https://ja.wikipedia.org/wiki/ガンマ分布) $Ga(a)$, $Ga(k,\theta)$

```Python:Python3
# ガンマ分布 Ga(a)
a = 5
plt.figure()
r = gamma.rvs(a=a, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 20, 0.1)
y = gamma.pdf(x, a=a)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = gamma.stats(a=a, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# ガンマ分布 Ga(k,theta)
k = 5
theta = 1
plt.figure()
r = gamma.rvs(a=k, scale=theta, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 20, 0.1)
y = gamma.pdf(x, a=k, scale=theta)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = gamma.stats(a=a, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_gamma.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/b3a56bcb-f17a-787b-12a1-a9a7f3cd2247.png)


- [指数分布](https://ja.wikipedia.org/wiki/指数分布) $Exp(\lambda)$

```Python:Python3
# 指数分布 Exp(L)
L = 1
plt.figure()
r = expon.rvs(loc=0, scale=1/L, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 5, 0.1)
y = expon.pdf(x, loc=0, scale=1/L)   # pdf = Lexp(-Lx)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = expon.stats(loc=0, scale=1/L, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_exp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/36cbc638-a2cc-dc80-bf14-ee907393b7e7.png)


- [対数正規分布](https://ja.wikipedia.org/wiki/対数正規分布) $LN(\mu,\sigma)$

```Python:Python3
# 対数正規分布 LN(mu,sigma)
mu = 0
sigma = 1
plt.figure()
r = lognorm.rvs(s=sigma, scale=np.exp(mu), size=1000, random_state=0)
plt.hist(r, bins=100, density=True, label='histgram')
x = np.arange(0, 10, 0.05)
y = lognorm.pdf(x, s=sigma, scale=np.exp(mu))
plt.plot(x, y, label='pdf')
plt.xlim(0,10)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = lognorm.stats(s=sigma, scale=np.exp(mu), moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

logr = np.log(r)
plt.figure()
plt.hist(logr, bins=50, density=True, label='Log(LN)')
plt.legend()
x = np.arange(-5, 5, 0.1)
y = norm.pdf(x, loc=mu, scale=sigma)
plt.plot(x, y, label='N(0,1)_pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(logr):{:.2f}, var(logr):{:.2f}, std(logr):{:.2f}'.format(np.mean(logr),np.var(logr),np.std(logr)))
```

![histgram_lognorm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/3fc63c6b-26be-1252-8a1b-c129c0b60d97.png) ![histgram_log_lognorm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/e2a6f655-c03b-34ac-22a9-505b3e169678.png)


**R**

```R:R
# 一様分布 U(a,b)
a <- 10
b <- 20
set.seed(0)
r <- runif(n=1000, min=a, max=b)
hist(r, breaks=50)

# 正規分布 N(mu,sigma)
m = 0
s = 1
set.seed(0)
r <- rnorm(n=1000, mean=m, sd=s)
hist(r, breaks=50)

# t分布 t(n)
n = 10
set.seed(0)
r <- rt(n=1000, df=n)
hist(r, breaks=50)

# カイ二乗分布 X2(n)
n = 10
set.seed(0)
r <- rchisq(n=1000, df=n)
hist(r, breaks=50)

# F分布 F(n,m)
n = 5
m = 20
set.seed(0)
r <- rf(n=1000, df1=n, df2=m)
hist(r, breaks=50)

# ベータ分布 Be(a,b)
a = 2
b = 5
set.seed(0)
r <- rbeta(n=1000, shape1=a, shape2=b)
hist(r, breaks=50)

# ガンマ分布 Ga(k,theta)
k = 5
theta = 1
set.seed(0)
r <- rgamma(n=1000, shape=k, scale=theta)
hist(r, breaks=50)

# 指数分布 Exp(L)
L = 1
set.seed(0)
r <- rexp(n=1000, rate=L)
hist(r, breaks=50)

# 対数正規分布 LN(mu,sigma)
mu = 0
sigma = 1
set.seed(0)
r <- rlnorm(n=1000, meanlog=m, sdlog=s)
hist(r, breaks=50)

hist(log(r), breaks=50)
```

**VBA**

```VBA:VBA
Dim size As Integer
Dim i As Integer
Dim r() As Single
Dim a As Double
Dim b As Double
Dim m As Double
Dim s As Double
Dim n As Integer
Dim L As Double

size = 1000
ReDim r(1 To size)

' 一様乱数 U(a,b)
a = 10
b = 20
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = (b - a) * Rnd + a
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 2).value = r(i)
Next i

' 正規分布 N(0,1)
m = 0
s = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Norm_Inv(Rnd, m, s)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 3).value = r(i)
Next i

't分布 t(10)
n = 10
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.T_Inv(Rnd, n)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 4).value = r(i)
Next i

' カイ二乗分布 Chi2(10)
n = 10
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.ChiSq_Inv(Rnd, n)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 5).value = r(i)
Next i

' F分布 F(5,20)
n = 5
m = 20
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.F_Inv(Rnd, n, m)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 6).value = r(i)
Next i

' ベータ分布 Be(2,5)
a = 2
b = 5
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Beta_Inv(Rnd, a, b)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 7).value = r(i)
Next i

' ガンマ分布 Ga(5,1)
a = 5
b = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Gamma_Inv(Rnd, a, b)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 8).value = r(i)
Next i

' 指数 Exp(1)=Ga(1,1)
L = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Gamma_Inv(Rnd, 1, L)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 9).value = r(i)
Next i
'なし

' 対数正規分布 LN(0,1)
m = 5
s = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.LogNorm_Inv(Rnd, m, s)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 10).value = r(i)
Next i
```


## 離散確率分布に従う乱数
次に、[離散確率分布](https://ja.wikipedia.org/wiki/離散確率分布)（continuous probability distribution）に従う乱数の発生を説明します。

### 離散確率分布に従う乱数の発生

**Python**

```Python:Python3
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom

n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
plt.figure()
plt.hist(r)
plt.show()
```
![histgram_binom1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/14152d2f-808b-5a81-16e3-d25f882460ec.png)

乱数の発生はこれでできていますが、ヒストグラムの描き方について補足します。
乱数は0～10の11種類発生していますが、このヒストグラムは10個に分かれています。
離散確率分布で整数しか出てきませんので、このようにビン（ヒストグラムの棒）の数を自動で設定すると変な感じになることがあります。そこで、ビンの数を直接指定します。また、x軸の目盛りも1刻みで表示されるように設定します。

```Python:Python3
n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
xbins = np.arange(0, n+1+1, 1)
print(xbins)
# [ 0  1  2  3  4  5  6  7  8  9 10 11]
plt.figure()
plt.hist(r, bins=xbins)
plt.xticks(xbins)
plt.show()
```
![histgram_binom2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/229202ab-897b-4313-7c06-56a192b61032.png)

これでもよいのですが、
- $x = 0$ の度数は、x軸の目盛り0の上ではなく、0と1の間に
- $x = 1$ の度数は、x軸の目盛り1の上ではなく、1と2の間に
- $x = 2$ の度数は、x軸の目盛り2の上ではなく、2と3の間に
・・・となっているのが見にくいので、グラフ全体を左にずらします。

```Python:Python3
n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
xbins = np.arange(0, n+1+1, 1)
print(xbins)
# [ 0  1  2  3  4  5  6  7  8  9 10 11]
plt.figure()
plt.hist(r, bins=xbins, align='left')
plt.xticks(xbins)
plt.show()
```
![histgram_binom3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/86591ecf-16d6-8327-fb87-4bcdc7bbb589.png)

最後に、確率質量関数（pmf）を重ねて描きます。確率質量関数は全体が1になるように正規化されていますので、ヒストグラムも全体を1になるように正規化しています。

```Python:Python3
n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
xbins = np.arange(0, n+1+1, 1)
print(xbins)
# [ 0  1  2  3  4  5  6  7  8  9 10 11]
plt.figure()
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
plt.xticks(xbins)
x = np.arange(0, n+1, 1)
y = binom.pmf(x, n, p)
plt.plot(x, y, marker='o', label='pdf')
plt.legend()
plt.show()
```
![histgram_binom4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/82cacae8-4ca7-63bf-5377-21b86db05623.png)


**R**

```R:R
# 二項分布
n <- 10
p <- 0.5
set.seed(0)
r <- rbinom(n=1000, size=n, prob=p)
hist(r)

# ビンを設定
hist(r, breaks=0:(n+1))

# ビンの区分はデフォルトでは right-closed (left open) intervals なので
# left-closed (right open) intervals に
hist(r, breaks=0:(n+1), right=F)

# 全体を1になるように正規化
hist(r, breaks=0:(n+1), right=F, prob=T)
x <- seq(0, 11, by=1)
y <- dbinom(x, size=n, prob=p)
par(new=T)
plot(x, y, xlab="", ylab="", col="red", type="l")
```

**VBA**
整数1から6の一様乱数

```VBA:VBA
' 一様乱数（整数）1, ... , 6
Rnd -1
Randomize 1
For i = 1 To 100
    Debug.Print Int(6 * Rnd + 1)
Next i
```


### 離散確率分布の例

**Python**

```Python:Python3
from scipy.stats import bernoulli, binom, poisson, nbinom, geom, hypergeom, randint
```

- [ベルヌーイ試行](https://ja.wikipedia.org/wiki/ベルヌーイ試行) $Bernoulli(p)$

```Python:Python3
# ベルヌーイ試行 Bernoulli(p)
p = 0.3
plt.figure()
xbins = np.arange(0, 2+1, 1)
r = bernoulli.rvs(p=p, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, 1+1, 1)
y = bernoulli.pmf(x, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = bernoulli.stats(p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_bernoulli.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/89b0cb2f-05ef-7abf-d991-fdadf62d4b6d.png)


- [二項分布](https://ja.wikipedia.org/wiki/二項分布) $B(n,p)$

```Python:Python3
# 二項分布 B(n,p)
# n回の独立なベルヌーイ試行Bernoulli(p)を行ったときの成功回数の分布
n = 10
p = 0.3
plt.figure()
xbins = np.arange(0, n+1+1, 1)
r = binom.rvs(n=n, p=p, size=1000, random_state=0)
plt.hist(r, bins=xbin, align='left', density=True, label='histgram')
x = np.arange(0, n+1, 1)
y = binom.pmf(x, n=n, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = binom.stats(n=n, p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_binom.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/bdbde8ba-7425-644f-5716-c54074cc1e90.png)


- [ポアソン分布](https://ja.wikipedia.org/wiki/ポアソン分布) $Po(\lambda)$

```Python:Python3
# ポアソン分布 Po(L)
# 所与の時間中に平均L回の頻度で発生する事象が発生する回数の分布
L = 5
plt.figure()
xbins = np.arange(0, L*3+1, 1)
r = poisson.rvs(mu=L, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, L*3+1, 1)
y = poisson.pmf(x, mu=L)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = poisson.stats(mu=L, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_poisson.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/a7be553e-fcd9-8e94-ccbe-05804645ceb4.png)


- [負の二項分布](https://ja.wikipedia.org/wiki/負の二項分布) $NB(n,p)$

```Python:Python3
# 負の二項分布(negative binomial distribution) NB(n,p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったときn回の成功を得るまでに失敗した試行回数の分布
n = 5
p = 0.4
plt.figure()
xbins = np.arange(0, n*5+1, 1)
r = nbinom.rvs(n=n, p=p, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, n*5+1, 1)
y = nbinom.pmf(x, n=n, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = nbinom.stats(n=n, p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_nbinom.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/e62cdfec-ed56-2ca7-0ecf-721318527750.png)


- [幾何分布](https://ja.wikipedia.org/wiki/幾何分布) $G(p)$

```Python:Python3
# 幾何分布(geometric distribution) G(p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったとき初めて成功するまでの試行回数の分布
# NB(1,p)+1
p = 0.4
plt.figure()
xbins = np.arange(1, 15+1, 1)
r = geom.rvs(p=p, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(1, 15+1, 1)
y = geom.pmf(x, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = geom.stats(p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_geom.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/02f785bc-245a-5c83-1d61-9bd284287852.png)


- [超幾何分布](https://ja.wikipedia.org/wiki/超幾何分布) $HG(M,n,N)$

```Python:Python3
# 超幾何分布(hypergeometric distribution) HG(M,n,N)
# n個の成功状態をもつM個の要素よりなる母集団からN個の要素を非復元抽出したときに含まれる成功状態の個数の分布
M = 20
n = 7
N = 12
plt.figure()
xbins = np.arange(0, min(n,N)+1, 1)
r = hypergeom.rvs(M=M, n=n, N=N, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, min(n,N)+1, 1)
y = hypergeom.pmf(x, M=M, n=n, N=N)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = hypergeom.stats(M=M, n=n, N=N, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_hypergeom.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/20d76953-2138-4409-1bdf-b902760cc4b2.png)


- 一様分布（整数）

```Python:Python3
# 一様分布（整数）
# a～bまでの整数の一様分布
a = 1
b = 6
plt.figure()
xbins = np.arange(a, b+1+1, 1)
r = randint.rvs(low=a, high=b+1, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(a, b+1, 1)
y = randint.pmf(x, low=a, high=b+1)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = randint.stats(low=a, high=b+1, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```

![histgram_randint.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/800718/a9b991bd-dc4c-9e61-ea0e-1ba841fe01e6.png)


**R**

```R:R
# ベルヌーイ試行 Bernoulli(p)
p <- 0.3
set.seed(0)
r <- rbinom(n=1000, size=1, prob=p)
hist(r, breaks=0:(1+1), right=F, prob=T)

# 二項分布 B(n,p)
# n回の独立なベルヌーイ試行Bernoulli(p)を行ったときの成功回数の分布
n <- 10
p <- 0.3
set.seed(0)
r <- rbinom(n=1000, size=n, prob=p)
hist(r, breaks=0:(n+1), right=F, prob=T)

# ポアソン分布 Po(L)
# 所与の時間中に平均L回の頻度で発生する事象が発生する回数の分布
L <- 5
set.seed(0)
r <- rpois(n=1000, lambda=L)
hist(r, breaks=0:20, right=F, prob=T)

# 負の二項分布(negative binomial distribution) NB(n,p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったときn回の成功を得るまでに失敗した試行回数の分布
n <- 5
p <- 0.4
set.seed(0)
r <- rnbinom(n=1000, size=n, prob=p)
hist(r, breaks=0:40, right=F, prob=T)

# 幾何分布(geometric distribution) G(p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったとき初めて成功するまでの試行回数の分布
# NB(1,p)+1
p <- 0.4
set.seed(0)
r <- rgeom(n=1000, prob=p)
hist(r, breaks=0:20, right=F, prob=T)

# 超幾何分布(hypergeometric distribution) HG(M,n,N)
# n個の成功状態をもつM個の要素よりなる母集団からN個の要素を非復元抽出したときに含まれる成功状態の個数の分布
M <- 20
n <- 7
N <- 12
set.seed(0)
r <- rhyper(nn=1000, n=M, m=n, k=N) # 乱数の個数はnn
hist(r, breaks=0:(min(n, N)+1), right=F, prob=T)

# 一様分布（整数）
# a～bまでの整数の一様分布
a <- 1
b <- 6
set.seed(0)
r <- sample(1:6, size=1000, replace=T)
r <- sample.int(n=6, size=1000, replace=T)
hist(r, breaks=a:(b+1), right=F, prob=T)
```

**VBA**
VBAでは、`Rnd`関数で一様乱数が発生できます。それ以外の分布に従う乱数は、[WorksheetFunction](https://docs.microsoft.com/ja-jp/office/vba/api/excel.worksheetfunction)の分布関数の逆関数を使えば発生できます。Excelの統計関数については[こちら](https://support.microsoft.com/ja-jp/office/統計関数-参照-624dac86-a375-4435-bc25-76d659719ffd)。

```VBA:VBA
Dim size As Integer
Dim i As Integer
Dim r() As Single
Dim a As Double
Dim b As Double
Dim m As Double
Dim s As Double
Dim n As Integer
Dim L As Double

size = 1000
ReDim r(1 To size)

' 一様乱数 U(a,b)
a = 10
b = 20
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = (b - a) * Rnd + a
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 2).value = r(i)
Next i

' 正規分布 N(0,1)
m = 0
s = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Norm_Inv(Rnd, m, s)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 3).value = r(i)
Next i

't分布 t(10)
n = 10
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.T_Inv(Rnd, n)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 4).value = r(i)
Next i

' カイ二乗分布 Chi2(10)
n = 10
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.ChiSq_Inv(Rnd, n)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 5).value = r(i)
Next i

' F分布 F(5,20)
n = 5
m = 20
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.F_Inv(Rnd, n, m)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 6).value = r(i)
Next i

' ベータ分布 Be(2,5)
a = 2
b = 5
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Beta_Inv(Rnd, a, b)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 7).value = r(i)
Next i

' ガンマ分布 Ga(5,1)
a = 5
b = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Gamma_Inv(Rnd, a, b)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 8).value = r(i)
Next i

' 指数 Exp(1)=Ga(1,1)
L = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Gamma_Inv(Rnd, 1, L)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 9).value = r(i)
Next i
'なし

' 対数正規分布 LN(0,1)
m = 5
s = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.LogNorm_Inv(Rnd, m, s)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 10).value = r(i)
Next i

' 一様乱数（整数）
a = 1
b = 6
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = Int((b - a + 1) * Rnd + a)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 18).value = r(i)
Next i
```


# まとめ
### 一覧
各言語で使用する関数等を一覧にまとめます。比較のために、EXCELでの計算も示しました。

***乱数***

||Python|R|VBA|EXCEL|
|:--|:--|:--|:--|:--|
|[一様分布](https://ja.wikipedia.org/wiki/一様分布) <br>$U(a,b)$|`uniform.rvs(loc=a, scale=b-a, size=N)`|`runif(n=N, min=a, max=b)`|`(b - a) * Rnd + a`|`=(b-a)*RAND()+b`|
|[正規分布](https://ja.wikipedia.org/wiki/正規分布) <br>$N(m,s^2)$|`norm.rvs(loc=m, scale=s, size=N)`|`rnorm(n=N, mean=m, sd=s)`|`WorksheetFunction.`<br>`Norm_Inv(Rnd, m, s)`|`=NORM.INV(`<br>`RAND(),m,s)`|
|[t分布](https://ja.wikipedia.org/wiki/T分布) <br>$t(n)$|`t.rvs(df=n, size=N)`|`rt(n=N, df=n)`|`WorksheetFunction.`<br>`T_Inv(Rnd, n)`|`=T.INV(`<br>`RAND(),n)`|
|[カイ二乗分布](https://ja.wikipedia.org/wiki/カイ二乗分布) <br>$\chi^2(n)$|`chi2.rvs(df=n, size=N)`|`rchisq(n=N, df=n)`|`WorksheetFunction.`<br>`ChiSq_Inv(Rnd, n)`|`=CHISQ.INV(`<br>`RAND(),n)`|
|[F分布](https://ja.wikipedia.org/wiki/F分布) <br>$F(n,m)$|`f.rvs(dfn=n, dfd=m, size=N)`|`rf(n=N, df1=n, df2=m)`|`WorksheetFunction.`<br>`F_Inv(Rnd, n, m)`|`=F.INV(`<br>`RAND(),n,m)`|
|[ベータ分布](https://ja.wikipedia.org/wiki/ベータ分布) <br>$Be(a,b)$|`beta.rvs(a=a, b=b, size=N)`|`rbeta(n=N, shape1=a, shape2=b)`|`WorksheetFunction.`<br>`Beta_Inv(Rnd, a, b)`|`=BETA.INV(`<br>`RAND(),a,b)`|
|[ガンマ分布](https://ja.wikipedia.org/wiki/ガンマ分布) <br>$Ga(k,theta)$|`gamma.rvs(a=k, scale=theta, size=N)`|`rgamma(n=N, shape=k, scale=theta)`|`WorksheetFunction.`<br>`Gamma_Inv(Rnd, a, b)`|`=GAMMA.INV(`<br>`RAND(),a,b)`|
|[指数分布](https://ja.wikipedia.org/wiki/指数分布) <br>$Exp(L)$|`expon.rvs(loc=0, scale=1/L, size=N)`|`rexp(n=N, rate=L)`|`WorksheetFunction.`<br>`Gamma_Inv(Rnd, 1, L)`|`=GAMMA.INV(`<br>`RAND(),1,L)`|
|[対数正規分布](https://ja.wikipedia.org/wiki/対数正規分布) <br>$LN(mu,s)$|`lognorm.rvs(s=s, scale=np.exp(m), size=N)`|`rlnorm(n=N, meanlog=m, sdlog=s)`|`WorksheetFunction.`<br>`LogNorm_Inv(Rnd, m, s)`|`=LOGNORM.INV(`<br>`RAND(),m,s)`|
|[ベルヌーイ試行](https://ja.wikipedia.org/wiki/ベルヌーイ試行) <br>$Bernoulli(p)$|`bernoulli.rvs(p=p, size=N)`|`rbinom(n=N, size=1, prob=p)`|||
|[二項分布](https://ja.wikipedia.org/wiki/二項分布) <br>$B(n,p)$|`binom.rvs(n=n, p=p, size=N)`|`rbinom(n=N, size=n, prob=p)`|||
|[ポアソン分布](https://ja.wikipedia.org/wiki/ポアソン分布) <br>$Po(L)$|`poisson.rvs(mu=L, size=N)`|`rpois(n=N, lambda=L)`||||
|[負の二項分布](https://ja.wikipedia.org/wiki/負の二項分布) <br>$NB(n,p)$|`nbinom.rvs(n=n, p=p, size=N)`|`rnbinom(n=N, size=n, prob=p)`|||
|[幾何分布](https://ja.wikipedia.org/wiki/幾何分布) <br>$G(p)$|`geom.rvs(p=p, size=N)`|`rgeom(n=N, prob=p)`|||
|[超幾何分布](https://ja.wikipedia.org/wiki/超幾何分布) $HG(M,n,N)$|`hypergeom.rvs(M=M, n=n, N=N, size=N)`|`rhyper(nn=N, n=M, m=n, k=N)`|||
|一様分布（整数）|`randint.rvs(low=a, high=b+1, size=N)`|`sample(1:6, size=N, replace=T)`<br>`sample.int(n=6, size=N, replace=T)`|`Int((b - a + 1) * Rnd + a)`|`=INT((b-a+1)*RAND()+a)`|
注意）
Pythonでは、`scipy.stats`から、必要な関数をインポートしたうえで使用。
`from scipy.stats import uniform, norm, chi2, t, f, beta, gamma, lognorm, expon`、`from scipy.stats import bernoulli, binom, poisson, nbinom, geom, hypergeom, randint`


**Pythonのライブラリの比較**

||Python(random)|Python(numpy.random)|Python(scipy.stats)|
|:--|:--|:--|:--|
|import|`import random`|`import numpy as np`|`from scipy.stats import norm`等|
|seedの固定|`random.seed(0)`|`np.random.seed(0)`|関数の引数で`random_state=0`|
|[一様分布](https://ja.wikipedia.org/wiki/一様分布) <br>$U(0,1)$|`random.random()`|`np.random.`<br>`rand(N)`<br>`np.random.`<br>`random_sample(size=N)`<br>`np.random.`<br>`random(size=N)`<br>`np.random.`<br>`ranf(size=N)`<br>`np.random.`<br>`sample(size=N)`|`uniform.rvs(loc=a, scale=b-a, size=N)`|
|[一様分布](https://ja.wikipedia.org/wiki/一様分布) <br>$U(a,b)$|`random.`<br>`uniform(a, b)`|`np.random.`<br>`uniform(low=a, high=b, size=N)`|`uniform.rvs(loc=a, scale=b-a, size=N)`|
|標準正規分布 <br>$N(0,1)$|`random.`<br>`gauss(0, 1)`<br>`random.`<br>`normalvariate(0, 1)`|`np.random.`<br>`randn(N)`<br>`np.random.`<br>`standard_normal(size=N)`|`norm.rvs(loc=m, scale=s, size=N)`|
|[正規分布](https://ja.wikipedia.org/wiki/正規分布) <br>$N(m,s^2)$|`random.`<br>`gauss(mu=m, sigma=s)`<br>`random.`<br>`normalvariate(mu=m, sigma=s)`|`np.random.`<br>`normal(loc=m, scale=s, size=N)`|`norm.rvs(loc=m, scale=s, size=N)`|
|[t分布](https://ja.wikipedia.org/wiki/T分布) <br>$t(n)$||`np.random.`<br>`standard_t(df=n, size=N)`|`t.rvs(df=n, size=N)`|
|[カイ二乗分布](https://ja.wikipedia.org/wiki/カイ二乗分布) <br>$\chi^2(n)$||`np.random.`<br>`chisquare(df=n, size=N)`|`chi2.rvs(df=n, size=N)`|
|[F分布](https://ja.wikipedia.org/wiki/F分布) <br>$F(n,m)$||`np.random.`<br>`f(dfnum=n, dfden=m, size=N)`|`f.rvs(dfn=n, dfd=m, size=N)`|
|[ベータ分布](https://ja.wikipedia.org/wiki/ベータ分布) <br>$Be(a,b)$|`random.`<br>`betavariate(alpha, beta)`|`np.random.`<br>`beta(a, b, size=N)`|`beta.rvs(a=a, b=b, size=N)`|
|[ガンマ分布](https://ja.wikipedia.org/wiki/ガンマ分布) <br>$Ga(k,theta)$|`random.`<br>`gammavariate(alpha, beta)`|`np.random.`<br>`gamma(shape, scale, size=N)`|`gamma.rvs(a=k, scale=theta, size=N)`|
|[指数分布](https://ja.wikipedia.org/wiki/指数分布) <br>$Exp(L)$|`random.`<br>`expovariate(lambd)`|`np.random.`<br>`exponential(scale=1/L, size=N)`|`expon.rvs(loc=0, scale=1/L, size=N)`|
|[対数正規分布](https://ja.wikipedia.org/wiki/対数正規分布) <br>$LN(mu,s)$|`random.`<br>`lognormvariate(mu=m, sigma=s)`|`np.random.`<br>`lognormal(mean=m, sigma=s, size=N)`|`lognorm.rvs(s=s, scale=np.exp(m), size=N)`|
|[ベルヌーイ試行](https://ja.wikipedia.org/wiki/ベルヌーイ試行) <br>$Bernoulli(p)$|||`bernoulli.rvs(p=p, size=N)`|
|[二項分布](https://ja.wikipedia.org/wiki/二項分布) <br>$B(n,p)$||`np.random.`<br>`binomial(n=n, p=p, size=N)`|`binom.rvs(n=n, p=p, size=N)`|
|[ポアソン分布](https://ja.wikipedia.org/wiki/ポアソン分布) <br>$Po(L)$||`np.random.`<br>`poisson(lam=L, size=N)`|`poisson.rvs(mu=L, size=N)`|
|[負の二項分布](https://ja.wikipedia.org/wiki/負の二項分布) <br>$NB(n,p)$||`np.random.`<br>`negative_binomial(n=n, p=p, size=N)`|`nbinom.rvs(n=n, p=p, size=N)`|
|[幾何分布](https://ja.wikipedia.org/wiki/幾何分布) <br>$G(p)$||`np.random.`<br>`geometric(p=p, size=N)`|`geom.rvs(p=p, size=N)`|
|[超幾何分布](https://ja.wikipedia.org/wiki/超幾何分布) $HG(M,n,N)$||`np.random.`<br>`hypergeometric(ngood, nbad, nsample, size=N)`|`hypergeom.rvs(M=M, n=n, N=N, size=N)`|
|一様分布（整数）|`random.`<br>`randint(a,b)`<br>`random.`<br>`randrange(start=a, stop=b, step=1)`|`np.random.`<br>`randint(low=a, high=b+1, size=N)`<br>`np.random.`<br>`random_integers(low=a, high=b, size=N)`|`randint.rvs(low=a, high=b+1, size=N)`|
|サンプリング（１個）|`random.`<br>`choice(seq=x)`|||
|サンプリング（重複あり）復元抽出|`random.`<br>`choices(population=x, k=N)`|`np.random.`<br>`choice(a=x, size=N, replace=True, p)`||
|サンプリング（重複なし）非復元抽出|`random.`<br>`sample(population=x, k=N)`|`np.random.`<br>`choice(a=x, size=N, replace=False, p)`||
|シャッフル|`random.`<br>`shuffle(x)`|`np.random.`<br>`shuffle(x)`<br>`np.random.`<br>`permutation(x)`|
注意）指数分布の`random.expovariate(lambd)`の引数`lambd`は $Exp(\lambda)$ の $\lambda$ （平均の値の逆数）で、本来は "lambda" だが、Pythonの予約語で使えないためこのような表記になっている。


### プログラム全体
参考までに使ったプログラムの全体を示します。

**Python**

```Python:Python3
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

m = 0
s = 1
r = norm.rvs(m, s, size=1000, random_state=0)
plt.figure()
plt.hist(r, bins=50)
plt.show()

m = 0
s = 1
r = norm.rvs(m, s, size=1000, random_state=0)
plt.figure()
plt.hist(r, bins=50, density=True, label='histgram') # 全体を1に正規化
x = np.arange(-5, 5, 0.1)
y = norm.pdf(x, m, s)
plt.plot(x, y, label='pdf')
plt.legend()
plt.show()

from scipy.stats import uniform, norm, chi2, t, f, beta, gamma, lognorm, expon

# 一様分布 U(a,b)
a = 10
b = 20
plt.figure()
r = uniform.rvs(loc=a, scale=b-a, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(a, b, 0.1)
y = uniform.pdf(x, loc=a, scale=b-a)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = uniform.stats(loc=a, scale=b-a, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 正規分布 N(m,s)
m = 0
s = 1
plt.figure()
r = norm.rvs(loc=m, scale=s, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(-5, 5, 0.1)
y = norm.pdf(x, loc=m, scale=s)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = norm.stats(loc=m, scale=s, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# t分布 t(n)
n = 10
plt.figure()
r = t.rvs(df=n, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(-5, 5, 0.1)
y = t.pdf(x, df=n)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = t.stats(df=n, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# カイ二乗分布 X2(n)
n = 10
plt.figure()
r = chi2.rvs(df=n, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 30, 0.1)
y = chi2.pdf(x, df=n)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = chi2.stats(df=n, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# F分布 F(n,m)
n = 5
m = 20
plt.figure()
r = f.rvs(dfn=n, dfd=m, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 5, 0.1)
y = f.pdf(x, dfn=n, dfd=m)   # numerator(分子), denominator(分母)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = f.stats(dfn=n, dfd=m, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# ベータ分布 Be(a,b)
a = 2
b = 5
plt.figure()
r = beta.rvs(a=a, b=b, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 1, 0.01)
y = beta.pdf(x, a=a, b=b)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = beta.stats(a=a, b=b, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# ガンマ分布 Ga(a)
a = 5
plt.figure()
r = gamma.rvs(a=a, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 20, 0.1)
y = gamma.pdf(x, a=a)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = gamma.stats(a=a, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# ガンマ分布 Ga(k,theta)
k = 5
theta = 1
plt.figure()
r = gamma.rvs(a=k, scale=theta, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 20, 0.1)
y = gamma.pdf(x, a=k, scale=theta)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = gamma.stats(a=a, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 指数分布 Exp(L)
L = 1
plt.figure()
r = expon.rvs(loc=0, scale=1/L, size=1000, random_state=0)
plt.hist(r, bins=50, density=True, label='histgram')
x = np.arange(0, 5, 0.1)
y = expon.pdf(x, loc=0, scale=1/L)   # pdf = Lexp(-Lx)
plt.plot(x, y, label='pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = expon.stats(loc=0, scale=1/L, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 対数正規分布 LN(mu,sigma)
mu = 0
sigma = 1
plt.figure()
r = lognorm.rvs(s=sigma, scale=np.exp(mu), size=1000, random_state=0)
plt.hist(r, bins=100, density=True, label='histgram')
x = np.arange(0, 10, 0.05)
y = lognorm.pdf(x, s=sigma, scale=np.exp(mu))
plt.plot(x, y, label='pdf')
plt.xlim(0,10)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = lognorm.stats(s=sigma, scale=np.exp(mu), moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

logr = np.log(r)
plt.figure()
plt.hist(logr, bins=50, density=True, label='Log(LN)')
plt.legend()
x = np.arange(-5, 5, 0.1)
y = norm.pdf(x, loc=mu, scale=sigma)
plt.plot(x, y, label='N(0,1)_pdf')
plt.legend(loc='upper right')
plt.show()
print('mean(logr):{:.2f}, var(logr):{:.2f}, std(logr):{:.2f}'.format(np.mean(logr),np.var(logr),np.std(logr)))


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom

n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
plt.figure()
plt.hist(r)
plt.show()

n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
xbins = np.arange(0, n+1+1, 1)
print(xbins)
# [ 0  1  2  3  4  5  6  7  8  9 10 11]
plt.figure()
plt.hist(r, bins=xbins)
plt.xticks(xbins)
plt.show()

n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
xbins = np.arange(0, n+1+1, 1)
print(xbins)
# [ 0  1  2  3  4  5  6  7  8  9 10 11]
plt.figure()
plt.hist(r, bins=xbins, align='left')
plt.xticks(xbins)
plt.show()

n = 10
p = 0.5
r = binom.rvs(n, p, size=1000, random_state=0)
xbins = np.arange(0, n+1+1, 1)
print(xbins)
# [ 0  1  2  3  4  5  6  7  8  9 10 11]
plt.figure()
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
plt.xticks(xbins)
x = np.arange(0, n+1, 1)
y = binom.pmf(x, n, p)
plt.plot(x, y, marker='o', label='pdf')
plt.legend()
plt.show()

from scipy.stats import bernoulli, binom, poisson, nbinom, geom, hypergeom, randint

# ベルヌーイ試行 Bernoulli(p)
p = 0.3
plt.figure()
xbins = np.arange(0, 2+1, 1)
r = bernoulli.rvs(p=p, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, 1+1, 1)
y = bernoulli.pmf(x, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = bernoulli.stats(p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 二項分布 B(n,p)
# n回の独立なベルヌーイ試行Bernoulli(p)を行ったときの成功回数の分布
n = 10
p = 0.3
plt.figure()
xbins = np.arange(0, n+1+1, 1)
r = binom.rvs(n=n, p=p, size=1000, random_state=0)
plt.hist(r, bins=xbin, align='left', density=True, label='histgram')
x = np.arange(0, n+1, 1)
y = binom.pmf(x, n=n, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = binom.stats(n=n, p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# ポアソン分布 Po(L)
# 所与の時間中に平均L回の頻度で発生する事象が発生する回数の分布
L = 5
plt.figure()
xbins = np.arange(0, L*3+1, 1)
r = poisson.rvs(mu=L, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, L*3+1, 1)
y = poisson.pmf(x, mu=L)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = poisson.stats(mu=L, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 負の二項分布(negative binomial distribution) NB(n,p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったときn回の成功を得るまでに失敗した試行回数の分布
n = 5
p = 0.4
plt.figure()
xbins = np.arange(0, n*5+1, 1)
r = nbinom.rvs(n=n, p=p, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, n*5+1, 1)
y = nbinom.pmf(x, n=n, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = nbinom.stats(n=n, p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 幾何分布(geometric distribution) G(p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったとき初めて成功するまでの試行回数の分布
# NB(1,p)+1
p = 0.4
plt.figure()
xbins = np.arange(1, 15+1, 1)
r = geom.rvs(p=p, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(1, 15+1, 1)
y = geom.pmf(x, p=p)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = geom.stats(p=p, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 超幾何分布(hypergeometric distribution) HG(M,n,N)
# n個の成功状態をもつM個の要素よりなる母集団からN個の要素を非復元抽出したときに含まれる成功状態の個数の分布
M = 20
n = 7
N = 12
plt.figure()
xbins = np.arange(0, min(n,N)+1, 1)
r = hypergeom.rvs(M=M, n=n, N=N, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(0, min(n,N)+1, 1)
y = hypergeom.pmf(x, M=M, n=n, N=N)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = hypergeom.stats(M=M, n=n, N=N, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))

# 一様分布（整数）
# a～bまでの整数の一様分布
a = 1
b = 6
plt.figure()
xbins = np.arange(a, b+1+1, 1)
r = randint.rvs(low=a, high=b+1, size=1000, random_state=0)
plt.hist(r, bins=xbins, align='left', density=True, label='histgram')
x = np.arange(a, b+1, 1)
y = randint.pmf(x, low=a, high=b+1)
plt.plot(x, y, marker='o', label='pmf')
plt.xticks(xbins)
plt.legend(loc='upper right')
plt.show()
print('mean(r):{:.2f}, var(r):{:.2f}, std(r):{:.2f}'.format(np.mean(r),np.var(r),np.std(r)))
mean, var, skew, kurt = randint.stats(low=a, high=b+1, moments='mvsk')
print('mean:{:.2f}, var:{:.2f}, skew:{:.2f}, kurt:{:.2f}'.format(mean,var,skew,kurt))
```


**R**

```R:R
m <- 0
s <- 1
set.seed(0) # 乱数のシードをセット
r <- rnorm(n=1000, mean=m, sd=s)
hist(r)
hist(r, breaks=50)

m <- 0
s <- 1
set.seed(0)
r <- rnorm(n=1000, mean=m, sd=s)
xrange <- c(-4,4)
yrange <- c(0,0.6)
hist(r, breaks=50, prob=T, xlim=xrange, ylim=yrange, ann=F)
x <- seq(-5, 5, by=0.1)
y <- dnorm(x, mean=0, sd=1)
par(new=T)
plot(x, y, xlim=xrange, ylim=yrange, xlab="", ylab="", col="red", type="l")
par(new=T)
curve(dnorm, from=-4, to=4, ylim=yrange, col="blue", ylab="")


# 一様分布 U(a,b)
a <- 10
b <- 20
set.seed(0)
r <- runif(n=1000, min=a, max=b)
hist(r, breaks=50)

# 正規分布 N(mu,sigma)
m = 0
s = 1
set.seed(0)
r <- rnorm(n=1000, mean=m, sd=s)
hist(r, breaks=50)

# t分布 t(n)
n = 10
set.seed(0)
r <- rt(n=1000, df=n)
hist(r, breaks=50)

# カイ二乗分布 X2(n)
n = 10
set.seed(0)
r <- rchisq(n=1000, df=n)
hist(r, breaks=50)

# F分布 F(n,m)
n = 5
m = 20
set.seed(0)
r <- rf(n=1000, df1=n, df2=m)
hist(r, breaks=50)

# ベータ分布 Be(a,b)
a = 2
b = 5
set.seed(0)
r <- rbeta(n=1000, shape1=a, shape2=b)
hist(r, breaks=50)

# ガンマ分布 Ga(k,theta)
k = 5
theta = 1
set.seed(0)
r <- rgamma(n=1000, shape=k, scale=theta)
hist(r, breaks=50)

# 指数分布 Exp(L)
L = 1
set.seed(0)
r <- rexp(n=1000, rate=L)
hist(r, breaks=50)

# 対数正規分布 LN(mu,sigma)
mu = 0
sigma = 1
set.seed(0)
r <- rlnorm(n=1000, meanlog=m, sdlog=s)
hist(r, breaks=50)

hist(log(r), breaks=50)


# 二項分布
n <- 10
p <- 0.5
set.seed(0)
r <- rbinom(n=1000, size=n, prob=p)
hist(r)

# ビンを設定
hist(r, breaks=0:(n+1))

# ビンの区分はデフォルトでは right-closed (left open) intervals なので
# left-closed (right open) intervals に
hist(r, breaks=0:(n+1), right=F)

# 全体を1になるように正規化
hist(r, breaks=0:(n+1), right=F, prob=T)
x <- seq(0, 11, by=1)
y <- dbinom(x, size=n, prob=p)
par(new=T)
plot(x, y, xlab="", ylab="", col="red", type="l")


# ベルヌーイ試行 Bernoulli(p)
p <- 0.3
set.seed(0)
r <- rbinom(n=1000, size=1, prob=p)
hist(r, breaks=0:(1+1), right=F, prob=T)

# 二項分布 B(n,p)
# n回の独立なベルヌーイ試行Bernoulli(p)を行ったときの成功回数の分布
n <- 10
p <- 0.3
set.seed(0)
r <- rbinom(n=1000, size=n, prob=p)
hist(r, breaks=0:(n+1), right=F, prob=T)

# ポアソン分布 Po(L)
# 所与の時間中に平均L回の頻度で発生する事象が発生する回数の分布
L <- 5
set.seed(0)
r <- rpois(n=1000, lambda=L)
hist(r, breaks=0:20, right=F, prob=T)

# 負の二項分布(negative binomial distribution) NB(n,p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったときn回の成功を得るまでに失敗した試行回数の分布
n <- 5
p <- 0.4
set.seed(0)
r <- rnbinom(n=1000, size=n, prob=p)
hist(r, breaks=0:40, right=F, prob=T)

# 幾何分布(geometric distribution) G(p)
# 独立なベルヌーイ試行Bernoulli(p)を行ったとき初めて成功するまでの試行回数の分布
# NB(1,p)+1
p <- 0.4
set.seed(0)
r <- rgeom(n=1000, prob=p)
hist(r, breaks=0:20, right=F, prob=T)

# 超幾何分布(hypergeometric distribution) HG(M,n,N)
# n個の成功状態をもつM個の要素よりなる母集団からN個の要素を非復元抽出したときに含まれる成功状態の個数の分布
M <- 20
n <- 7
N <- 12
set.seed(0)
r <- rhyper(nn=1000, n=M, m=n, k=N) # 乱数の個数はnn
hist(r, breaks=0:(min(n, N)+1), right=F, prob=T)

# 一様分布（整数）
# a～bまでの整数の一様分布
a <- 1
b <- 6
set.seed(0)
r <- sample(1:6, size=1000, replace=T)
r <- sample.int(n=6, size=1000, replace=T)
hist(r, breaks=a:(b+1), right=F, prob=T)
```


**VBA**

```VBA:VBA
Sub test_rand()

Randomize 1
For i = 1 To 10
    Debug.Print Rnd
Next i

Dim n As Integer
Dim r() As Single
n = 10
ReDim r(1 To n)
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = Rnd
    Debug.Print i, r(i)
Next i
' 1             0.3335753
' 2             6.816387E-02
' 3             0.5938293
' 4             0.7660395
' 5             0.1892894
' 6             0.5373986
' 7             0.3269944
' 8             0.393937
' 9             7.341915E-02
' 10            0.8315425

' 一様乱数（整数）1, ... , 6
Rnd -1
Randomize 1
For i = 1 To 10
    Debug.Print Int(6 * Rnd + 1)
Next i

End Sub

Sub test_rand_examples()
Dim size As Integer
Dim i As Integer
Dim r() As Single
Dim a As Double
Dim b As Double
Dim m As Double
Dim s As Double
Dim n As Integer
Dim L As Double

size = 1000
ReDim r(1 To size)

' 一様乱数 U(a,b)
a = 10
b = 20
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = (b - a) * Rnd + a
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 2).value = r(i)
Next i

' 正規分布 N(0,1)
m = 0
s = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Norm_Inv(Rnd, m, s)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 3).value = r(i)
Next i

't分布 t(10)
n = 10
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.T_Inv(Rnd, n)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 4).value = r(i)
Next i

' カイ二乗分布 Chi2(10)
n = 10
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.ChiSq_Inv(Rnd, n)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 5).value = r(i)
Next i

' F分布 F(5,20)
n = 5
m = 20
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.F_Inv(Rnd, n, m)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 6).value = r(i)
Next i

' ベータ分布 Be(2,5)
a = 2
b = 5
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Beta_Inv(Rnd, a, b)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 7).value = r(i)
Next i

' ガンマ分布 Ga(5,1)
a = 5
b = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Gamma_Inv(Rnd, a, b)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 8).value = r(i)
Next i

' 指数 Exp(1)=Ga(1,1)
L = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.Gamma_Inv(Rnd, 1, L)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 9).value = r(i)
Next i
'なし

' 対数正規分布 LN(0,1)
m = 5
s = 1
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = WorksheetFunction.LogNorm_Inv(Rnd, m, s)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 10).value = r(i)
Next i

' 一様乱数（整数）
a = 1
b = 6
Rnd -1
Randomize 1
For i = LBound(r) To UBound(r)
    r(i) = Int((b - a + 1) * Rnd + a)
    Debug.Print i, r(i)
    ActiveSheet.Cells(i + 1, 18).value = r(i)
Next i

End Sub
```


# 参考
- Python
   - [Python チュートリアル](https://docs.python.org/ja/3/tutorial/index.html)
   - [SciPy.org](https://scipy.org/)
      - [Numpy and Scipy Documentation](https://docs.scipy.org/doc/)
      - [SciPy Reference](https://docs.scipy.org/doc/scipy/reference/)
         - API Reference
             - [Statistical functions (scipy.stats)](https://docs.scipy.org/doc/scipy/reference/stats.html)
         - [SciPy Tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/index.html)
             - [Statistics (scipy.stats)](https://docs.scipy.org/doc/scipy/reference/tutorial/stats.html)
   - [random](https://docs.python.org/ja/3/library/random.html)
   - [numpy.random](https://docs.scipy.org/doc/numpy-1.12.0/reference/routines.random.html)
- R
   - [RjpWiki](http://www.okadajp.org/RWiki/)
   - [R-Tips](http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html)
      - [60. 確率分布と乱数](http://cse.naro.affrc.go.jp/takezawa/r-tips/r/60.html)
      - [61. データの分布とヒストグラム・密度推定](http://cse.naro.affrc.go.jp/takezawa/r-tips/r/61.html)
- VBA
   - [Excel VBA リファレンス](https://docs.microsoft.com/ja-jp/office/vba/api/overview/excel)
      - [WorksheetFunction オブジェクト](https://docs.microsoft.com/ja-jp/office/vba/api/excel.worksheetfunction)
         - [メソッド](https://docs.microsoft.com/ja-jp/office/vba/api/excel.worksheetfunction#methods)
   - [Visual Basic 言語リファレンス](https://docs.microsoft.com/ja-jp/office/vba/language/reference/user-interface-help/visual-basic-language-reference)
      - [関数 (Visual Basic for Applications)](https://docs.microsoft.com/ja-jp/office/vba/language/reference/functions-visual-basic-for-applications)
      - [タスク別のキーワード](https://docs.microsoft.com/ja-jp/office/vba/language/reference/user-interface-help/keywords-by-task)
      - [オブジェクト (Visual Basic for Applications)](https://docs.microsoft.com/ja-jp/office/vba/language/reference/objects-visual-basic-for-applications)
      - [ステートメント](https://docs.microsoft.com/ja-jp/office/vba/language/reference/statements)
- Excel 
   - リファレンス
      - [統計関数 (リファレンス)](https://support.microsoft.com/ja-jp/office/統計関数-参照-624dac86-a375-4435-bc25-76d659719ffd)