はじめに
機械学習の勉強を始めたVBAユーザです。
備忘録としてPython・Rの文法をVBAと比較しながらまとめています。
今回は、乱数についてまとめます。
確率分布については、別記事にしました。
目次
乱数
今回は、乱数についてまとめます。
連続確率分布に従う乱数
まず、連続確率分布（continuous probability distribution）に従う乱数の発生から説明します。
連続確率分布に従う乱数の発生
例として、正規分布に従う乱数で説明します。
Python
Pythonでは、科学技術計算ライブラリの SciPy の統計関係のサブパッケージ scipy.stats を使います。scipy.stats から使用する関数をインポートしておきます。（なお、標準ライブラリのrandomモジュールや、Numpyのnumpy.randomモジュールを使っても乱数を発生できます。これについては下参照。）
計算に使うので、Numpy もインポートしておきます。また、ヒストグラムを描画するために Matplotlib の matplotlib.pyplot もインポートしておきます。
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()
乱数の発生はこれでできていますが、ヒストグラムの描き方について補足します。
plt.hist(r)でベクトル r のヒストグラムが描画できます。bins=50 でビンの数をしています。
これに、正規分布の確率密度関数（pdf）を重ねて描画してみます。確率密度関数は全体が1になるように正規化されていますので、ヒストグラムも全体が1になるように正規化しておきます。plt.histの引数でdensity=Trueと指定することで、全体を1に正規化したヒストグラムにできます。
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()
なお、Pythonでは他にも、Python標準ライブラリのrandomモジュールや、Numpyのnumpy.randomモジュールを使っても乱数を発生できます。
標準ライブラリのrandomモジュールでは1個の乱数を発生します。
Numpyのnumpy.randomモジュールでは、数ベクトルの形だけでなく配列の形の乱数も発生できます。
詳しくは下の一覧表をご覧ください。
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()
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
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
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を使った方が「一貫性のある文法で合理的に」描けます。（こちらの記事参照。）
VBA
VBAでは、0以上1未満の一様乱数ならRnd 関数で発生できます。引数の値によって乱数の生成方法が変わります（こちら）。
Randomize ステートメントで、Rnd 関数の乱数ジェネレーターを初期化します。引数に数値を指定するとその数値を新しいシード値として使用され、引数を省略するとシステムタイマーから返される値が新しいシード値として使用されます。
正規乱数など、他の分布に従う乱数については、下の例を参照。
VBA
Randomize 1
For i = 1 To 1000
Debug.Print Rnd
Next i
数値の引数を指定した Randomize を使用する直前に、負の引数を指定した Rnd を呼び出しておくことで、同じ乱数のシーケンスを発生させることができます。
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
Python3
from scipy.stats import uniform, norm, chi2, t, f, beta, gamma, lognorm, expon
- 一様分布 $U(a,b)$
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))
- 正規分布 $N(\mu,\sigma^2)$
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))
- t分布 $t(n)$
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))
- カイ二乗分布 $\chi^2(n)$
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))
- F分布 $F(n,m)$
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))
- ベータ分布 $Be(a,b)$
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))
- ガンマ分布 $Ga(a)$, $Ga(k,\theta)$
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))
- 指数分布 $Exp(\lambda)$
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))
- 対数正規分布 $LN(\mu,\sigma)$
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)))
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
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
離散確率分布に従う乱数
次に、離散確率分布（continuous probability distribution）に従う乱数の発生を説明します。
離散確率分布に従う乱数の発生
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()
乱数の発生はこれでできていますが、ヒストグラムの描き方について補足します。
乱数は0～10の11種類発生していますが、このヒストグラムは10個に分かれています。
離散確率分布で整数しか出てきませんので、このようにビン（ヒストグラムの棒）の数を自動で設定すると変な感じになることがあります。そこで、ビンの数を直接指定します。また、x軸の目盛りも1刻みで表示されるように設定します。
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()
これでもよいのですが、
- $x = 0$ の度数は、x軸の目盛り0の上ではなく、0と1の間に
- $x = 1$ の度数は、x軸の目盛り1の上ではなく、1と2の間に
- $x = 2$ の度数は、x軸の目盛り2の上ではなく、2と3の間に
・・・となっているのが見にくいので、グラフ全体を左にずらします。
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()
最後に、確率質量関数（pmf）を重ねて描きます。確率質量関数は全体が1になるように正規化されていますので、ヒストグラムも全体を1になるように正規化しています。
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()
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
' 一様乱数（整数）1, ... , 6
Rnd -1
Randomize 1
For i = 1 To 100
Debug.Print Int(6 * Rnd + 1)
Next i
離散確率分布の例
Python
Python3
from scipy.stats import bernoulli, binom, poisson, nbinom, geom, hypergeom, randint
- ベルヌーイ試行 $Bernoulli(p)$
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))
- 二項分布 $B(n,p)$
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))
- ポアソン分布 $Po(\lambda)$
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))
- 負の二項分布 $NB(n,p)$
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))
- 幾何分布 $G(p)$
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))
- 超幾何分布 $HG(M,n,N)$
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))
- 一様分布（整数）
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))
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の分布関数の逆関数を使えば発生できます。Excelの統計関数についてはこちら。
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 | |
|---|---|---|---|---|
|
一様分布 $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 |
|
正規分布 $N(m,s^2)$ |
norm.rvs(loc=m, scale=s, size=N) |
rnorm(n=N, mean=m, sd=s) |
WorksheetFunction.Norm_Inv(Rnd, m, s)
|
=NORM.INV(RAND(),m,s)
|
|
t分布 $t(n)$ |
t.rvs(df=n, size=N) |
rt(n=N, df=n) |
WorksheetFunction.T_Inv(Rnd, n)
|
=T.INV(RAND(),n)
|
|
カイ二乗分布 $\chi^2(n)$ |
chi2.rvs(df=n, size=N) |
rchisq(n=N, df=n) |
WorksheetFunction.ChiSq_Inv(Rnd, n)
|
=CHISQ.INV(RAND(),n)
|
|
F分布 $F(n,m)$ |
f.rvs(dfn=n, dfd=m, size=N) |
rf(n=N, df1=n, df2=m) |
WorksheetFunction.F_Inv(Rnd, n, m)
|
=F.INV(RAND(),n,m)
|
|
ベータ分布 $Be(a,b)$ |
beta.rvs(a=a, b=b, size=N) |
rbeta(n=N, shape1=a, shape2=b) |
WorksheetFunction.Beta_Inv(Rnd, a, b)
|
=BETA.INV(RAND(),a,b)
|
|
ガンマ分布 $Ga(k,theta)$ |
gamma.rvs(a=k, scale=theta, size=N) |
rgamma(n=N, shape=k, scale=theta) |
WorksheetFunction.Gamma_Inv(Rnd, a, b)
|
=GAMMA.INV(RAND(),a,b)
|
|
指数分布 $Exp(L)$ |
expon.rvs(loc=0, scale=1/L, size=N) |
rexp(n=N, rate=L) |
WorksheetFunction.Gamma_Inv(Rnd, 1, L)
|
=GAMMA.INV(RAND(),1,L)
|
|
対数正規分布 $LN(mu,s)$ |
lognorm.rvs(s=s, scale=np.exp(m), size=N) |
rlnorm(n=N, meanlog=m, sdlog=s) |
WorksheetFunction.LogNorm_Inv(Rnd, m, s)
|
=LOGNORM.INV(RAND(),m,s)
|
|
ベルヌーイ試行 $Bernoulli(p)$ |
bernoulli.rvs(p=p, size=N) |
rbinom(n=N, size=1, prob=p) |
||
|
二項分布 $B(n,p)$ |
binom.rvs(n=n, p=p, size=N) |
rbinom(n=N, size=n, prob=p) |
||
|
ポアソン分布 $Po(L)$ |
poisson.rvs(mu=L, size=N) |
rpois(n=N, lambda=L) |
||
|
負の二項分布 $NB(n,p)$ |
nbinom.rvs(n=n, p=p, size=N) |
rnbinom(n=N, size=n, prob=p) |
||
|
幾何分布 $G(p)$ |
geom.rvs(p=p, size=N) |
rgeom(n=N, prob=p) |
||
| 超幾何分布 $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)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
|
|
一様分布 $U(0,1)$ |
random.random() |
np.random.rand(N)np.random.random_sample(size=N)np.random.random(size=N)np.random.ranf(size=N)np.random.sample(size=N)
|
uniform.rvs(loc=a, scale=b-a, size=N) |
|
一様分布 $U(a,b)$ |
random.uniform(a, b)
|
np.random.uniform(low=a, high=b, size=N)
|
uniform.rvs(loc=a, scale=b-a, size=N) |
| 標準正規分布 $N(0,1)$ |
random.gauss(0, 1)random.normalvariate(0, 1)
|
np.random.randn(N)np.random.standard_normal(size=N)
|
norm.rvs(loc=m, scale=s, size=N) |
|
正規分布 $N(m,s^2)$ |
random.gauss(mu=m, sigma=s)random.normalvariate(mu=m, sigma=s)
|
np.random.normal(loc=m, scale=s, size=N)
|
norm.rvs(loc=m, scale=s, size=N) |
|
t分布 $t(n)$ |
np.random.standard_t(df=n, size=N)
|
t.rvs(df=n, size=N) |
|
|
カイ二乗分布 $\chi^2(n)$ |
np.random.chisquare(df=n, size=N)
|
chi2.rvs(df=n, size=N) |
|
|
F分布 $F(n,m)$ |
np.random.f(dfnum=n, dfden=m, size=N)
|
f.rvs(dfn=n, dfd=m, size=N) |
|
|
ベータ分布 $Be(a,b)$ |
random.betavariate(alpha, beta)
|
np.random.beta(a, b, size=N)
|
beta.rvs(a=a, b=b, size=N) |
|
ガンマ分布 $Ga(k,theta)$ |
random.gammavariate(alpha, beta)
|
np.random.gamma(shape, scale, size=N)
|
gamma.rvs(a=k, scale=theta, size=N) |
|
指数分布 $Exp(L)$ |
random.expovariate(lambd)
|
np.random.exponential(scale=1/L, size=N)
|
expon.rvs(loc=0, scale=1/L, size=N) |
|
対数正規分布 $LN(mu,s)$ |
random.lognormvariate(mu=m, sigma=s)
|
np.random.lognormal(mean=m, sigma=s, size=N)
|
lognorm.rvs(s=s, scale=np.exp(m), size=N) |
|
ベルヌーイ試行 $Bernoulli(p)$ |
bernoulli.rvs(p=p, size=N) |
||
|
二項分布 $B(n,p)$ |
np.random.binomial(n=n, p=p, size=N)
|
binom.rvs(n=n, p=p, size=N) |
|
|
ポアソン分布 $Po(L)$ |
np.random.poisson(lam=L, size=N)
|
poisson.rvs(mu=L, size=N) |
|
|
負の二項分布 $NB(n,p)$ |
np.random.negative_binomial(n=n, p=p, size=N)
|
nbinom.rvs(n=n, p=p, size=N) |
|
|
幾何分布 $G(p)$ |
np.random.geometric(p=p, size=N)
|
geom.rvs(p=p, size=N) |
|
| 超幾何分布 $HG(M,n,N)$ |
np.random.hypergeometric(ngood, nbad, nsample, size=N)
|
hypergeom.rvs(M=M, n=n, N=N, size=N) |
|
| 一様分布（整数） |
random.randint(a,b)random.randrange(start=a, stop=b, step=1)
|
np.random.randint(low=a, high=b+1, size=N)np.random.random_integers(low=a, high=b, size=N)
|
randint.rvs(low=a, high=b+1, size=N) |
| サンプリング（１個） |
random.choice(seq=x)
|
||
| サンプリング（重複あり）復元抽出 |
random.choices(population=x, k=N)
|
np.random.choice(a=x, size=N, replace=True, p)
|
|
| サンプリング（重複なし）非復元抽出 |
random.sample(population=x, k=N)
|
np.random.choice(a=x, size=N, replace=False, p)
|
|
| シャッフル |
random.shuffle(x)
|
np.random.shuffle(x)np.random.permutation(x)
|
|
注意）指数分布のrandom.expovariate(lambd)の引数lambdは $Exp(\lambda)$ の $\lambda$ （平均の値の逆数）で、本来は "lambda" だが、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
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
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
- R
- VBA
- Excel
- リファレンス