Python で対数正規分布の計算を実行する方法を調べてみました。
結果を記載します。


対数正規分布とは

統計学入門 - 東京大学出版会 6.11 対数正規分布には以下記載されています。

ランダムに世帯を選びその年間所得Xを調べると、低い方は一定限度があるが、高い方には明確な限度がない。このような場合、logX を考えるのが自然である(図6.15)。そこで、logX が正規分布に従うならば、もとのXは対数正規分布 log-normal distribution に従うという。

現実で、対数正規分布と近い分布をとる例には以下があります。

  • オプション: プレミアムの評価方法
  • 低所得者範囲での所得分布

Python での対数正規分布の計算方法

対数正規分布はlog-normal distribution呼ばれます。log-normal distribution python検索した結果以下が見つかりました。

numpy.random.lognormal を使う

numpy.random.lognormal対数正規分布に従う配列を生成できます。

import numpy as np
# mu: 基礎となる正規分布の平均値、 sigma: 基礎となる正規分布の標準偏差    
mu, sigma = 3., 1. 
s = np.random.lognormal(mu, sigma, 1000)
print(s)

[ 53.65039347  16.42483065   5.21176677   4.73970538   5.4776618
   8.21022889  34.58747085 121.00853537  13.86098751   6.43636124
   4.06565487  31.76621297   9.4816259    7.10916311  33.48658321
   2.05610029  50.77300074  28.96997044  13.51811062  98.36436465
  .....
  44.15313676  38.45845665   4.91379163  21.96371365   7.26879496
  18.76733173  14.52764593  38.79054122  10.20387406  50.45289727
  51.39174082  54.2082318  205.19325454   5.61990921   2.62172464]

ヒストグラムを描画します。

import matplotlib.pyplot as plt
import numpy as np
# mu: 基礎となる正規分布の平均値、 sigma: 基礎となる正規分布の標準偏差    
mu, sigma = 3., 1. 
s = np.random.lognormal(mu, sigma, 1000)
count, bins, ignored = plt.hist(s, 100, density=True, align='mid')
plt.show()

201904080311_output_4_0.png - Google ドライブ

mu sigma変更すると、グラフの分布が変わります。

import matplotlib.pyplot as plt
import numpy as np
# mu: 基礎となる正規分布の平均値、 sigma: 基礎となる正規分布の標準偏差    
mu, sigma = 6., 1.5 
s = np.random.lognormal(mu, sigma, 1000)
count, bins, ignored = plt.hist(s, 100, density=True, align='mid')
plt.show()

201904080311_output_6_0.png - Google ドライブ

scipy.stats.lognorm を使う

scipy.stats.lognorm使うと、対数正規分布に従う配列の他、確率密度関数等の生成ができます。

from scipy.stats import lognorm
s = 0.954
# モーメントの計算をする
mean, var, skew, kurt = lognorm.stats(s, moments='mvsk')
print(mean, var, skew, kurt)

1.576264803741382 3.6886797556399684 5.464256148333118 81.30583448651502

import matplotlib.pyplot as plt
import numpy as np
# shpae (形状パラメータ)に 0.954 を指定する   
s = 0.954
# %点関数 1%から 99% 点までを100等分した配列を生成する   
x = np.linspace(lognorm.ppf(0.01, s),lognorm.ppf(0.99, s), 100)
fig, ax = plt.subplots(1, 1)
# pdf (Probability density function) 確率密度関数
ax.plot(x, lognorm.pdf(x, s),'r-', lw=5, alpha=0.6, label='lognorm pdf')
# shapeを指定して、変数で受ける    
rv = lognorm(s)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
r = lognorm.rvs(s, size=1000)
# ヒストグラムに描画する    
ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()

201904080311_output_9_0.png - Google ドライブ

import matplotlib.pyplot as plt
import numpy as np
# shpae (形状パラメータ)に 0.5を指定する    
s = 0.5
# %点関数 1%から 99% 点までを100等分した配列を生成する   
x = np.linspace(lognorm.ppf(0.01, s),lognorm.ppf(0.99, s), 100)
fig, ax = plt.subplots(1, 1)
# pdf (Probability density function) 確率密度関数
ax.plot(x, lognorm.pdf(x, s),'r-', lw=5, alpha=0.6, label='lognorm pdf')
# shapeを指定して、変数で受ける    
rv = lognorm(s)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
r = lognorm.rvs(s, size=1000)
# ヒストグラムに描画する    
ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
ax.legend(loc='best', frameon=False)
plt.show()

201904080311_output_10_0.png - Google ドライブ

参考

以下、記事作成中に参考にしました。

以上です。

コメント