データ分析系プログラマーのブログ

主にPythonを使ったデータ分析や機械学習をやっています。

連続一様分布について

確率密度関数が常に一定の値をとる分布のことを一様分布と言います。一様分布には、離散型確率分布と連続型確率分布があります。この記事では、連続一様分布について紹介しています。

連続一様分とは

確率変数がどのような値でも、その時の確率密度関数が一定の値をとる分布のことを連続一様分布といいます。

連続一様分布に従う確率密度関数と期待値と分散

確率変数Xがa <= X <= bにおける連続一様分布に従うとき、確率密度関数と期待値と分散はそれぞれ次のようになります。


f(x) = \frac{1}{b-a}\quad (a \leq X \leq b)

f(x) = 0\quad (X < a, X > b)

E(X) = \frac{a+b}{2}

V(X) = \frac{(b-a)^2}{12}

例えば、確率変数Xが1 <= X <= 4における連続一様分布に従うとき、確率密度関数と期待値と分散はそれぞれ以下のようになります。


f(x) = \frac{1}{b-a}\quad (1 \leq X \leq 4) = \frac{1}{4-1} =\frac{1}{3}

E(X) = \frac{a+b}{2}=\frac{1+4}{2}=\frac{5}{2}

V(X) = \frac{(b-a)^2}{12}=\frac{(4-1)^2}{12}=\frac{9}{12}=\frac{3}{4}

連続一様分布に従う累積分布関数

確率密度関数f(x)が連続一様分布に従う場合の累積分布関数F(x)は以下のようになります。


F(x) = \frac{x-a}{b-a}\quad (a \leq X < b)

F(x) = 0\quad (X < a)

F(x) = 1\quad (X \geq b)

例えば、2以上12以下の範囲で乱数を10,000個作成したとき、その中で6以下、もしくは9以上の値をとる確率は以下のようになります。


F(x) =(2 \leq X \leq 6)+(9 \leq X \leq 12)= \frac{6-2}{12-2}+\frac{12-9}{12-2}=\frac{4}{10}+\frac{3}{10}=\frac{7}{10}

連続一様分布に従う確率と期待値と分散の求め方その1

pythonで連続一様分布に従う確率、期待値、分散を求めるには、先ほどの式を使用して以下のように求めることができます。

from fractions import Fraction

a = 1
b = 4
fx = Fraction(1, (b - a))

ex = Fraction((a + b), 2)
vx = Fraction((b - a)**2, 12)

print("f(X ) = ", fx)
print("E(X) = ", ex)
print("V(X) = ", vx)

f(X ) = 1/3
E(X) = 5/2
V(X) = 3/4

連続一様分布に従う累積分布関数の求め方その1

pythonで連続一様分布に従う累積分布関数の確率を求める場合は、先ほどの式を使用して以下のように求めることができます。

a = 2
b = 12
x1 = 6
x2 = 9

print("F(X) = (2 < x < 6) + (9 < x < 12)なので")
Fx = Fraction((x1 - a), (b - a)) + Fraction((b - x2), (b - a))

print("F(X ) = ", Fx)

F(X) = (2 < x < 6) + (9 < x < 12)なので
F(X ) = 7/10

連続一様分布に従う期待値と分散の求め方その2(random.uniform)

random.uniformを使用すると、連続一様分布に従ったデータを生成することができます。以下の例では、random.uniformを使用して1から4までの連続一様分布に従うデータを10000個生成して、そのデータを使用して期待値と分散を求めています。

import random
from statistics import mean
from statistics import pvariance

x = []
for i in range(10000):
    x.append(random.uniform(1, 4))

m = mean(x)
v = pvariance(x)

print("E(X) = ", m)
print("V(X) = ", v)

E(X) = 2.5047683204489855
V(X) = 0.7423550594411364

連続一様分布に従う累積分布関数の求め方その2(random.uniform)

連続一様分布に従う累積分布関数の確率を求める場合は、以下のように条件に合うデータの数をカウントすることで確率を求めることができます。

x = []
for i in range(10000):
    x.append(random.uniform(2, 12))

p1 = 0
p2 = 0

for j in x:
    if j <= 6:
        p1 += 1 #6以下のデータをカウントする

for k in x:
    if k >= 9:
        p2 += 1 #9以上のデータをカウントする

p = (p1+ p2)/10000
print("F(X ) = ", p)

F(X ) = 0.6981

連続一様分布に従う期待値と分散の求め方その3(np.random.uniform)

np.random.uniformを使用すると、連続一様分布に従ったデータを生成することができます。以下の例では、np.random.uniformを使用して1から4までの連続一様分布に従うデータを10000個生成して、そのデータを使用して期待値と分散を求めています。

import numpy as np
from statistics import mean
from statistics import pvariance

x = np.random.uniform(1,4,10000)

m = mean(x)
v = pvariance(x)

print("E(X) = ", m)
print("V(X) = ", v)

E(X) = 2.50598766012386
V(X) = 0.7446555437307987

連続一様分布に従う累積分布関数の求め方その3(np.random.uniform)

連続一様分布に従う累積分布関数の確率を求める場合は、以下のように条件に合うデータの数をカウントすることで確率を求めることができます。

x = np.random.uniform(2,12,10000)

p1 = 0
p2 = 0

for i in x:
    if i <= 6:
        p1 += 1 #6以下のデーたをカウントする

for j in x:
    if j >= 9:
        p2 += 1 #9以上のデータをカウントする

p = (p1+ p2)/10000
print("F(X ) = ", p)

F(X ) = 0.702

連続一様分布に従う確率と期待値と分散の求め方その4(uniform)

scipy.statsのuniformを使用することでも、連続一様分布に従ったデータを生成することができます。以下の例では、scipy.statsのuniformを使用して1から4までの連続一様分布に従うデータを10000個生成して、そのデータを使用して期待値と分散を求めています。scipy.statsのuniformの場合、区間を指定するときは、loc, loc+scaleとなるように指定します。

from scipy.stats import uniform
from statistics import mean
from statistics import pvariance

x = np.linspace(uniform.ppf(0.01, loc=1.0, scale=3.0), uniform.ppf(0.99, loc=1.0, scale=3.0), 1)
p = uniform.pdf(x,loc=1, scale=3)

x2 = uniform.rvs(loc=1.0, scale=3.0, size=10000) #区間 loc, loc+scale

m = mean(x2)
v = pvariance(x2)

print("P(X) = ", p)
print("E(X) = ", m)
print("V(X) = ", v)

P(X) = [0.33333333]
E(X) = 2.4900863096103154
V(X) = 0.7406137576339912

連続一様分布に従う累積分布関数の求め方その4(uniform)

連続一様分布に従う累積分布関数の確率を求める場合は、以下のように条件に合うデータの数をカウントすることで確率を求めることができます。

x2 = uniform.rvs(loc=2.0, scale=10.0, size=10000) #区間 loc, loc+scale
p1 = 0
p2 = 0

for i in x2:
    if i <= 6:
        p1 += 1 #6以下のデータをカウントする

for j in x2:
    if j >= 9:
        p2 += 1 #9以上のデータをカウントする

p = (p1+ p2)/10000
print("F(X ) = ", p)

F(X ) = 0.7041

参考

統計学の時間 | 統計WEB

離散一様分布について

確率密度関数が常に一定の値をとる分布のことを一様分布と言います。一様分布には、離散型確率分布と連続型確率分布があります。この記事では、離散一様分布について紹介しています。

離散一様分布とは

離散一様分とは、サイコロのように離散型確率分布であり、1から6までの目が出る確率がすべて1/6というように同じ確率の分布を言います。

離散一様分布に従う確率と期待値と分散

確率変数Xが離散一様分布に従う時、Xの確率と期待値と分散は以下のようになります。


P(X = k) = \frac{1}{N}\quad (k = 1,2,3,\dots,N)

E(X) = \frac{N+1}{2}

V(X) = \frac{N^2+1}{12}

例えば、サイコロの出る目の確率は離散一様分布に従いますが、サイコロの出る目をXとすると、確率、期待値、分散は以下のようになります。


P(X = 1)=P(X = 2)=P(X = 3)=P(X = 4)=P(X = 5)=P(X = 6)=\frac{1}{6}

E(X)=\frac{6+1}{2}=3.5

V(X)=\frac{6^2-1}{12}=2.9

離散一様分布に従う確率と期待値と分散の求め方その1

pythonで離散一様分布に従う確率、期待値、分散を求めるには、先ほどの式を使用して以下のように求めることができます。

n = 6
px = 1/n
ex = (n + 1)/2
vx = (n**2 - 1)/12

print("P(X = 6) = ", px)
print("E(X) = ", ex)
print("V(X) = ", vx)

P(X = 6) = 0.16666666666666666
E(X) = 3.5
V(X) = 2.9166666666666665

離散一様分布に従う確率と期待値と分散の求め方その2

pythonで確率、期待値、分散を求める場合には、以下のように1から6までのデータを用意して求めることもできます。

import math

dice = [1,2,3,4,5,6]#サイコロの事象
ex = 0#期待値を演算するための変数
p = 1/len(dice)#サイコロの確率
sum01 = 0

for x in dice:
    px = (p*x)
    ex += px

for i in dice:
    r = (i - ex)**2
    sum01 += r

vx =p*sum01

print("P(X = 6) = ", p)
print("E(X) = ", ex)
print("V(X) = ", vx)

P(X = 6) = 0.16666666666666666
E(X) = 3.5
V(X) = 2.9166666666666665

離散一様分布に従う確率と期待値と分散の求め方その3

データが用意してある場合の確率、期待値、分散の求め方については、statisticsのライブラリを使用する場合もあります。meanやpvarianceを使用することで期待値や分散をあらかじめ要したデータから簡単に求めることができます。

from statistics import mean
from statistics import pvariance
dice = [1,2,3,4,5,6]

s = [i for i in dice if i == 1]
p = len(s)/len(dice)
m = mean(dice)
v = pvariance(dice)

print("P(X = 6) = ", p)
print("E(X) = ", m)
print("V(X) = ", v)

P(X = 6) = 0.16666666666666666
E(X) = 3.5
V(X) = 2.9166666666666665

離散一様分布に従う確率と期待値と分散の求め方その4

以下の例は、random.randintを使用して1から6までの整数をランダムに10000個生成したデータを使用して、それぞれ確率、期待値、分散を求めています。結果を見てみると、概ね解析的に求めた数値を同じ結果になっていることがわかります。

import random
discrete_uniform = [random.randint(1, 6) for i in range(10000)]

s = [i for i in discrete_uniform if i == 1]
p = len(s)/len(discrete_uniform)
m = mean(discrete_uniform)
v = pvariance(discrete_uniform)

print("P(X = 6) = ", p)
print("E(X) = ", m)
print("V(X) = ", v)

P(X = 6) = 0.1683
E(X) = 3.4941
V(X) = 2.92116519

参考

統計学の時間 | 統計WEB

指数分布について

連続型確率分布の一つに、機械が故障してから次に故障するまでの期間や、災害が起こってから次に起こるまでの期間のように、次に何かが起こるまでの期間が従う分布の指数分布があります。この記事では、pythonを使った指数分布の求め方について紹介しています。

指数分布とは

ある期間に平均して$\lambda$(ラムダ)回起こる現象が、次に起こるまでの期間$X$が指数分布に従うとき、$X=x$となる確率密度関数$f(x)$は次の式で表されます。


f(x) = \left\{
\begin{array}{ll}
\lambda e^{-\lambda x} & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.

例えば、1時間に平均20人のお客さんが送るお店に、あるお客がきて次のお客が来るまでの時間が10分となる確率密度は以下のように計算します。この場合は、10分は$\frac{1}{6}$時間となります。


f(\frac{1}{6})=\lambda e^{-\lambda x} = 20\times e^{-20\times \frac{1}{6}} = 0.714

指数分布に従う期待値と分散

確率変数$X$が指数分布に従う時、$X$の期待値と分散は以下のようになります。


E(X) = \frac{1}{\lambda}

V(X) = \frac{1}{\lambda^2}

例えば、先ほどの例の場合、$\lambda$は20となるので、期待値と分散は以下になります。


E(X) = \frac{1}{\lambda} =\frac{1}{20}

V(X) = \frac{1}{\lambda^2} =\frac{1}{400}

指数分布に従う累積分布関数

積分布関数は、ある期間に平均して$\lambda$回起こる現象が次に起こるまでの期間を$X$としたとき、期間$X$が$x$以下となる確率のことを言い、次のようになります。


F(x) = P (X \leq x) = 1 - e^{-\lambda x}

例えば、1時間に平均20人が来るお店に、ある客が来てから次の客が来るまでの時間が10分以内である確率は以下のように計算します。


F(\frac{1}{6}) = P (X \leq \frac{1}{6}) = 1 - e^{-20\times \frac{1}{6}} = 1 - 0.03569 = 0.96431

指数分布に従う確率密度の求め方

Pythonを使って指数分布に従う確率密度を求める場合は、以下のようになります。

import math

e = math.e
lambda1 = 20
x = 1/6

fx = lambda1*e**(-lambda1*x)
print("f(1/6) = ", fx)

f(1/6) = 0.7134798669450483

指数分布に従う累積分布関数の求め方その1

Pythonを使って、指数分布に従う累積分布関数の確率を求める場合は、以下のように解析的に求める方法があります。

e = math.e
lambda1 = 20
x = 1/6

Fx = 1 - e**(-lambda1*x)
print("P(X <= 1/6) = ", Fx)

P(X <= 1/6) = 0.9643260066527476

指数分布に従う累積分布関数の求め方その2(np.random.exponential)

Numpyのnp.random.exponentialを使用すると、指数分布に従うデータを生成することができます。なので、以下のように10分以内(1/6)の確率は以下のように計算することができます。

import numpy as np
size = 100000
lambda1 = 20
x = 1/6
beta = 1/lambda1

Fx = sum(np.random.exponential(beta, size) <= x)/size
print("P(X <= 1/6) = ", Fx)

P(X <= 1/6) = 0.96435

指数分布に従う累積分布関数の求め方その3(expon.rvs)

scipyのexpon.rvsを使用すると、指数分布に従うデータを生成することができます。なので、以下のように10分以内(1/6)の確率は以下のように計算することができます。

from scipy.stats import expon
size = 100000
lambda1 = 20
x = 1/6
beta = 1/lambda1

Fx = sum(expon.rvs(scale=beta, size=size) <= x)/size
print("P(X <= 1/6) = ", Fx)

P(X <= 1/6) = 0.96468

参考

統計学の時間 | 統計WEB

標準化について

正規分布に従う確率変数は、標準化を行うことによって標準正規分布に従う確率変数に変換することができます。標準化したデータは、異なる分布のデータ同士を比較したりすることができます。この記事では、pythonを使った標準化の方法について紹介しています。

標準化とは

ある確率変数Xが平均mu、分散sigma2正規分布に従う時、から平均muを引いて標準偏差sigmaで割った値をzとおくと、zは「平均が0、分散が1の標準正規分布」に従います。このような計算を標準化と呼び、以下のようにあらわします。標準化を行うことで単位や平均値が異なるデータ同士を比較できるようになります。


z = \frac{X-\mu}{\sigma}

例えば、あるクラスの数学と国語のテストの結果が以下のような場合、 ・数学 平均点:55点 標準偏差:15点 ・国語 平均点:45点 標準偏差:10点 A君は、数学が85点、国語が80点だった時、どちらの教科の方が順位が上かを比較するには以下のようにします。


z_m = \frac{X-\mu}{\sigma} = \frac{85-55}{15} = 2

z_j = \frac{X-\mu}{\sigma} = \frac{80-45}{10} = 3.5

この場合は、標準化した数値が大きい方が順位が上なので、国語の方が順位が上であることがわかります。

偏差値とは

偏差値は「平均が50点、標準偏差が10点」となるように、標準化した値zに10をかけて50を足したものです。


50 + 10 \times z

例えば、先ほどのA君の数学と国語の得点から偏差値を計算すると以下のようになります。

・数学の偏差値:50+10×2=70 ・国語の偏差値:50+10×3.5=85

Pythonを使った標準化の例その1(meanとpstdev)

以下の例では、random.randintを使って0から100までの整数をランダムに50個生成したデータを使って標準化の変換を行っています。標準化を行うためには、そのデータの平均と標準偏差を求める必要があります。以下では、statisticsのmeanとpstdevを使って平均と標準偏差を計算しています。このデータを標準化した結果、12番目にある88のデータは標準化することによって、約1.273となりました。また、このデータを点数とした場合、88の偏差値は約62.73と計算されました。

import random
data = [random.randint(0, 100) for i in range(50)]
print(data)
[28, 11, 73, 17, 44, 15, 4, 11, 60, 21, 100, 88, 28, 66, 61, 34, 16, 18, 66, 20, 95, 96, 100, 94, 76, 3, 81, 17, 30, 53, 84, 23, 97, 13, 89, 75, 50, 88, 16, 42, 40, 22, 11, 45, 26, 77, 42, 20, 86, 56]
from statistics import mean
from statistics import pstdev

mu = mean(data)
s = pstdev(data)
print("平均 = ", mu, "標準偏差", s)

平均 = 48.56 標準偏差 30.977514425789554

z = (data[11] - mu)/s
print(z)

1.2731815554305808

dev = 50 + 10*z
print(dev)

62.73181555430581

Pythonを使った標準化の例その2(zscore)

標準化の計算には、その他にもあります。例えば、scipyのzscoreを使用することによって、与えられたデータ(numpy形式の必要がある)から自動的に平均と標準偏差を求めて、標準化を行います。実際に標準化されたデータを確認してみると、88だったデータは、約1.273となりました。

from scipy.stats import zscore
arr_data = np.array(data)
data_std = zscore(arr_data)
print(data_std[11])

1.2731815554305808

Pythonを使った標準化の例その3(StandardScaler)

sklearnのStandardScalerを使うことでも、データを標準化することができます。StandardScalerは、機械学習を行う際の前処理などでもよく使用するライブラリです。この場合も、リスト形式のデータをnumpyに変換することで、StandardScalerを実行することができます。結果は他と同じように約1.273となっています。

import numpy as np
from sklearn.preprocessing import StandardScaler

arr_data = np.array(data)
arr_data = arr_data.reshape(-1, 1)
sc = StandardScaler()
data_std = sc.fit_transform(arr_data)
print(data_std[11])

[1.27318156]

標準化したデータを描画する

以下の例は、numpyによって生成したデータを使って、平均が50、標準偏差が10に従う正規分布の確率変数を標準化して描画しています。標準化の方法は、標準化の式を使用して、平均を50、標準偏差が10という設定で計算した場合のものがあります。その次の例では、StandardScalerを使って、標準化を行っていますが、この場合は、StandardScalerは与えられたデータから平均と標準偏差を求めて標準化するため、xは最初の例とは異なる値になってしまいます。先後の例は、scipyのnormを使って、平均50、標準偏差10の正規分布を生成し、xの値を標準化したデータで描画しています。この場合は、最初のグラフと同じ結果になります。

import numpy as np
import math
import matplotlib.pyplot as plt

x = np.linspace(0, 100, 10000)
e = math.e
pi = math.pi
mu = 50
s = 10

fx = (1/math.sqrt(2*pi)*s)*e**(-(((x - mu)**2)/(2*s**2)))
x_std = (x - mu)/s

fig, axes = plt.subplots(figsize = (12, 8))
plt.plot(x_std, fx)
plt.show()

f:id:hira03:20200211132335p:plain
statistics_standardization_distribution_e

import numpy as np
import math
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler

x = np.linspace(0, 100, 10000)
e = math.e
pi = math.pi
mu = 50
s = 10

fx = (1/math.sqrt(2*pi)*s)*e**(-(((x - mu)**2)/(2*s**2)))

x = x.reshape(-1, 1)
sc = StandardScaler()
x_std = sc.fit_transform(x)

fig, axes = plt.subplots(figsize = (12, 8))
plt.plot(x_std, fx)
plt.show()

f:id:hira03:20200211132411p:plain
statistics_StandardScaler_distribution_e

from scipy.stats import norm

x = np.linspace(0, 100, 10000)
mu = 50
s = 10

x_std = (x - mu)/s
fx = norm.pdf(x, loc=mu, scale=s)

fig, axes = plt.subplots(figsize = (12, 8))
plt.plot(x_std, fx*100)
plt.show()

f:id:hira03:20200211132442p:plain
statistics_standard_normal_distribution_e

参考

統計学の時間 | 統計WEB

正規分布について

連続型確率分布には正規分布と呼ばれる分布があります。正規分布は、検定や推定など様々な場面で用いられます。この記事では、pythonを使った正規分布ついて解析的に求める場合と、数値的に求める場合について紹介しています。

正規分布とは

正規分布とは、連続型確率分布の一つで、統計学における検定や推定、モデル作成などで用いられます。また、多くの統計的手法において、正規分布に従うことを仮定します。

正規分布に従う確率変数の確率密度関数の式

正規分布に従う確率変数Xの確率密度関数f(x)は次の式で表されます。


f(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}} \quad (-\infty < x < \infty)

または、eのかわりにexpを使うと以下のようになります。


f(x) = \frac{1}{\sqrt{2\pi}\sigma}exp\Bigl(-\frac{(x-\mu)^2}{2\sigma^2}\Bigr) \quad (-\infty < x < \infty)

標準正規分布に従う確率変数の確率密度関数の式

正規分布の中で、特に平均mu=0、分散sigma2=1である正規分布を「標準正規分布」といいます。先ほどの正規分布の式に平均mu=0、分散sigma2=1を代入すると以下のようになります。


f(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \quad (-\infty < x < \infty)

または


f(x) = \frac{1}{\sqrt{2\pi}}exp\Bigl(-\frac{x^2}{2}\Bigr) \quad (-\infty < x < \infty)

ネイピア数を使った正規分布のグラフの描き方

ここでは、ネイピア数を用いた正規分布を描画しています。やり方としては、numpyで0から100までの間に、10000の分割データを用意し、そのデータをネイピア数で表した正規分布の式で計算した結果をプロットしています。この時の正規分布は平均50、標準偏差10で計算しています。

import numpy as np
import math
import matplotlib.pyplot as plt

x = np.linspace(0, 100, 10000)
e = math.e
pi = math.pi
mu = 50
s = 10

fx = (1/math.sqrt(2*pi)*s)*e**(-(((x - mu)**2)/(2*s**2)))

fig, axes = plt.subplots(figsize = (12, 8))
plt.plot(x, fx)
plt.show()

f:id:hira03:20200210152606p:plain
statistics_normal_distribution_e

expを使った正規分布のグラフの描き方

expを用いた正規分布の場合は、ネイピア数の時と同じように、numpyでデータを用意します。expの場合は、numpyの一つ一つのデータに対して計算する必要があるので、for文で処理しています。

x = np.linspace(0, 100, 10000)
pi = math.pi
mu = 50
s = 10

x2 = -(((x - mu)**2)/(2*s**2))
fx = np.array([(1/math.sqrt(2*pi)*s)*math.exp(i) for i in x2])

fig, axes = plt.subplots(figsize = (12, 8))
plt.plot(x, fx)
plt.show()

f:id:hira03:20200210152638p:plain
statistics_normal_distribution_exp

random.gaussを使った正規分布の描き方

乱数を生成できるrandomにはガウス分布正規分布と同じ意味)を生成できるrandom.gaussがあります。これを使うと、正規分布にしたがったデータを生成することができます。random.gaussには平均と標準偏差の引数を設定できます。生成したデータをヒストグラムで可視化すると、先ほどの解析的に生成した正規分布のグラフとほぼ重なることがわかります。

import random
x = np.linspace(0, 100, 10000)
e = math.e
pi = math.pi
mu = 50
s = 10

fx = (1/math.sqrt(2*pi)*s)*e**(-(((x - mu)**2)/(2*s**2)))

normals = [random.gauss(mu, s) for i in range(10000)]

fig, axes = plt.subplots(figsize = (12, 8))
plt.hist(normals, 40, density = True)
plt.plot(x, 0.01*fx)
plt.show()

f:id:hira03:20200210152715p:plain
statistics_normal_distribution_gauss

random.normalvariateを使った正規分布の描き方

random.normalvariateを使うことでも、正規分布に従ったデータを生成することができます。生成したデータをヒストグラムで可視化すると、こちらも解析的に生成した正規分布のグラフとほぼ重なることがわかります。

import random
x = np.linspace(0, 100, 10000)
e = math.e
pi = math.pi
mu = 50
s = 10

fx = (1/math.sqrt(2*pi)*s)*e**(-(((x - mu)**2)/(2*s**2)))

normals = [random.normalvariate(mu, s) for i in range(10000)]

fig, axes = plt.subplots(figsize = (12, 8))
plt.hist(normals, 40, density = True)
plt.plot(x, 0.01*fx)
plt.show()

f:id:hira03:20200210152755p:plain
statistics_normal_distribution_normalvariate

numpy.random.normalを使った正規分布の描き方

numpyにも、乱数を生成する機能があります。numpy.random.normalを使うと、正規分布に従ったデータを生成することができます。numpy.random.normalの場合は、平均と標準偏差と生成するデータ数を指定するだけで正規分布に従うデータを生成することができます。

x = np.linspace(0, 100, 10000)
e = math.e
pi = math.pi
mu = 50
s = 10
size = 10000

fx = (1/math.sqrt(2*pi)*s)*e**(-(((x - mu)**2)/(2*s**2)))

normals = np.random.normal(mu , s, size)

fig, axes = plt.subplots(figsize = (12, 8))
plt.hist(normals, 40, density = True)
plt.plot(x, 0.01*fx)
plt.show()

f:id:hira03:20200210152840p:plain
statistics_normal_distribution_normal

scipy.stats.normを使った正規分布の描き方

scipy.stats.normを使用すると正規分布を簡単に描くことができます。この場合は、numpyで0から100までのデータを10000作り、norm.pdfで引数を設定することで正規分布の結果を得ることができます。必要な引数はnumpyデータと平均(loc)、標準偏差(scale)となっています。

from scipy.stats import norm

x = np.linspace(0, 100, 10000)
mu = 50
s = 10

fx = norm.pdf(x, loc=mu, scale=s)

fig, axes = plt.subplots(figsize = (12, 8))
plt.plot(x, fx)
plt.show()

f:id:hira03:20200210152909p:plain
statistics_normal_distribution_norm

参考

統計学の時間 | 統計WEB

Pythonを使った幾何分布の求め方

離散型確率分布には、成功確率がpの時のベルヌーイ試行を繰り返すとき、初めて成功するまでの試行回数に従う確率分布として幾何分布があります。この記事では、pythonを使った幾何分布の求め方について解析的に求める場合と、数値的に求める場合について紹介しています。

幾何分布とは

成功確率がpである独立なベルヌーイ試行を繰り返す時、初めて成功するまでの試行回数Xが従う確率分布を幾何分布といいます。


P(X = k) = (1-p)^{k-1} p

例えば、サイコロを2投目で初めて1が出る確率は以下のようになります。


P(X = 2) = (1-\frac{1}{6})^{2-1} \times\frac{1}{6} = \frac{5}{36} = 0.1388

幾何分布の期待値と分散

確率変数Xが成功確率pの幾何分布に従う時、Xの期待値と分散は以下のようになります。


E(X) = \frac{1}{p}\\
V(X) = \frac{1-p}{p^2}

例えば、けん玉を10回やった結果、4回は失敗する場合、次に失敗するまでの回数の期待値と分散は以下のようになります。


E(X) = \frac{1}{\frac{4}{10}} =\frac{5}{2} = 2.5\\
V(X) = \frac{1-\frac{4}{10}}{(\frac{4}{10})^2} = \frac{\frac{3}{5}}{\frac{4}{25}} = \frac{3}{5}\times\frac{25}{4} = \frac{75}{20} = 3.75

幾何分布を解析的に求める場合

幾何分布の確率を解析的に求める場合は、以下のように求めることができます。

p = 1/6
k = 2

px = p*(1 - p)**(k - 1)
print("P(X = 2) = ", px)

P(X = 2) = 0.1388888888888889

np.random.geometricを使って幾何分布を試行し数値的に確率を求める場合

np.random.geometricを使用すると、幾何分布を実際に試行することができます。以下の例では、確率が1/6の時の幾何分布を10000セット実行し、サイコロ2投目で1が出る回数を出力しています。そして、実際に実行したデータを使って幾何分布の確率を出しています。10000セットくらい実行したデータから求めた確率だと解析的に求めた確率に近くなっています。

import numpy as np
p = 1/6
k = 2

px = np.random.geometric(p,10000)
print("P(X = 2) = ", (px == k).sum() / 10000)

P(X = 2) = 0.1339

scipyのgeomを使って幾何分布を試行し数値的に確率を求める場合

scipyのgeomを使用することでも、幾何分布を実際に試行することができます。以下の例では、確率が1/6の時の幾何分布を10000セット実行し、サイコロ2投目で1が出る回数を出力しています。

from scipy.stats import geom
p = 1/6
k = 2

px = geom.rvs(p, size=10000)
print("P(X = 2) = ", (px == k).sum() / 10000)

P(X = 2) = 0.1383

Pythonを使ったポアソン分布の求め方

離散型確率分布のうち、結果が2つしかないような試行のことをベルヌーイ試行と言いますが、このうちnが大きく、pが小さい時をポアソン分布と言います。この記事では、pythonのライブラリを使ってポアソン分布の確率を解析的に求める場合と、数値的に求める場合について紹介しています。

ポアソン分布とは

ベルヌーイ試行のうちnが大きく、pが小さい時、「np = 一定(λ)」と考えることができ、ある期間にX回起きる確率分布をポアソン分布呼びます。


\lambda = np\\
e = 2.7182818 \dots\\
P(X = k) = \frac{e^{-\lambda}\lambda^k}{k!}

例えば、100個に1個の当たりのあるお菓子を10個買った時に、当たりのお菓子が2個含まれる場合の確率は以下のようになります。


\lambda = 10\times\frac{1}{100} = 0.1\\
e = 2.7182818 \dots\\
P(X = 2) = \frac{e^{-0.1}\times0.1^2}{2!} = 0.0045

ポアソン分布の期待値と分散

確率変数Xがポアソン分布に従う時、Xの期待値と分散は以下のようにλなります。


E(X) = \lambda\\
V(X) = \lambda

例えば、100個に1個の当たりのあるお菓子を10個買った時の期待値と分散はともに以下のようになります。


E(X) = \lambda = np = 10\times\frac{1}{100} = 0.1\\
V(X) = \lambda = 0.1

ポアソン分布を解析的に求める場合その1

ポアソン分布の確率を解析的に求める場合は、mathにネイピア数が用意されているので以下のように求めることができます。

import math

e = math.e
n = 10
p = 1/100
lambda1 = n*p
k = 2

px = (e**(-lambda1)*lambda1**k)/math.factorial(k)
print("P(X = 2) = ", px)

P(X = 2) = 0.004524187090179798

ポアソン分布を解析的に求める場合その2(poisson)

pythonのライブラリのscipyのpoissonを使用すると引数を指定するだけでポアソン分布の確率を求めることができます。

from scipy.stats import poisson

n = 10
p = 1/100
lambda1 = n*p
k = 2

px = poisson.pmf(k,lambda1)
print("P(X = 2) = ", px)

P(X = 2) = 0.004524187090179801

np.random.poissonを使ってポアソン分布を試行し数値的に確率を求める場合

np.random.poissonを使用すると、ポアソン分布を実際に試行することができます。以下の例では、ラムダが0.1の時のポアソン分布を500セット実行し、あたりの出る回数を出力しています。そして、実際に実行したデータを使ってポアソン分布の確率を出しています。100000セットくらい実行したデータから求めた確率だと解析的に求めた確率に近くなっています。

import numpy as np

n = 10
p = 1/100
lambda1 = n*p

s = np.random.poisson(lambda1, 500)

print(s)
[0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 0 0
 0 0 0 0 0 2 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 1
 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 1 0 1 0 0 0 0 2 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0
 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0
 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0]
s = sum(np.random.poisson(lambda1, 100000) == 2)/100000
print("P(X = 2) = ", s)

P(X = 2) = 0.00466

参考

統計学の時間 | 統計WEB