分布の不明な母集団からの1標本についてのメモ

はじめに

30代男性の毎月の飲み会代を調べることを考えたいと思います.
例えば23人にヒアリングするとします.
そうすることで、標本サイズ23の標本が1つ手に入ります.
ただし何らかの事情により標本は1度しか得られません.
本稿ではそうした状況についてのメモです.

母集団の分布について

今回のケースでは、観測者にとって母集団の分布が不明です.
という認識のもと、便宜上母集団の分布を与えます.
Bingにチャットしたところ、ガンマ分布なる分布が良いとのことでした.
個人的にガンマ分布のことをあまりよく知りません、そのためpythonコードとともに、グラフを以下に示します.

import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib #日本語化matplotlib
import seaborn as sns
sns.set(font="IPAexGothic") #日本語フォント設定
from scipy.stats import gamma


x = np.linspace(0, 100000, 1000)
y = gamma.pdf(x, a=2.5, scale=10000)

plt.plot(x, y)
plt.title('30代の毎月の飲み会代の分布')
plt.xlabel('飲み会代')
plt.ylabel('確率密度関数')
plt.show()

Pythonでガンマ分布を作成するには、SciPyライブラリのgamma関数を使用することができます。以下のコードは、30代の毎月の飲み会代の分布を作成する例です。この例では、aは形状パラメータで、scaleは尺度パラメータです。これらのパラメータを調整することで、ガンマ分布を調整することができます(Bing)

パラメータが2つあり、形状パラメータ、尺度パラメータという名前が存在するようです.
さて、非常に素晴らしい分布を手に入れることができました. 見通しをよくするために少し修正しました.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib #日本語化matplotlib
import seaborn as sns
sns.set(font="IPAexGothic") #日本語フォント設定
from scipy.stats import gamma

x = np.linspace(0, 100000, 1000)
y = gamma.pdf(x, a=2.5, scale=10000)

plt.plot(x, y)
plt.title('30代の毎月の飲み会代の分布')
plt.xlabel('飲み会代')
plt.ylabel('確率密度関数')

# 2.5%点と97.5%点を求める
low = gamma.ppf(0.025, a=2.5, scale=10000)
high = gamma.ppf(0.975, a=2.5, scale=10000)

# 信頼区間を描画する
plt.axvline(low, color='r', linestyle='--', label='2.5%点')
plt.axvline(high, color='r', linestyle='--', label='97.5%点')

plt.legend()

# グラフ内に数値を表示する
plt.text(low, 0, f'2.5%点: {low:.2f}', fontsize=10)
plt.text(high, 0, f'97.5%点: {high:.2f}', fontsize=10)

plt.show()

パッと見て理解しやすいようになりました.
Bingチャットのおかげで気が散らずに考えを進めることができ感謝です.

メモ

次に、この分布から標本サイズ23の標本を1度だけ得たいと思います.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib #日本語化matplotlib
import seaborn as sns
sns.set(font="IPAexGothic") #日本語フォント設定
from scipy.stats import gamma

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8), sharex=True)

x = np.linspace(0, 100000, 1000)
y = gamma.pdf(x, a=2.5, scale=10000)

ax1.plot(x, y, label='確率密度関数')
ax1.set_title('30代の毎月の飲み会代の分布')
ax1.set_xlabel('飲み会代')
ax1.set_ylabel('確率密度関数')

# 2.5%点と97.5%点を求める
low = gamma.ppf(0.025, a=2.5, scale=10000)
high = gamma.ppf(0.975, a=2.5, scale=10000)

# 信頼区間を描画する
ax1.axvline(low, color='r', linestyle='--', label='2.5%点')
ax1.axvline(high, color='r', linestyle='--', label='97.5%点')

# 標本をプロットする
np.random.seed(0) # 乱数生成器の状態を固定するために追加した行です。
sample = gamma.rvs(a=2.5, scale=10000, size=23)
ax1.scatter(sample, np.zeros_like(sample), color='black', label='標本')

# グラフ内に数値を表示する
ax1.text(low + 5000, gamma.pdf(low,a=2.5,scale=10000), f'2.5%点: {low:.2f}', fontsize=10)
ax1.text(high + 5000, gamma.pdf(high,a=2.5,scale=10000), f'97.5%点: {high:.2f}', fontsize=10)

sns.histplot(sample, ax=ax2)

# 標本平均を縦棒で描画する
mean = sample.mean()
ax2.axvline(mean, color='r', linestyle='--', label=f'標本平均: {mean:.2f}')
ax2.text(mean + 5000, gamma.pdf(high,a=2.5,scale=10000), f'標本平均: {mean:.2f}', fontsize=10)

# 凡例を表示する
ax1.legend()
plt.show()

上図の上側は先ほどのガンマ分布に標本をプロットしたものです.
上図の下側は標本をヒストグラムにしたもので、x軸は同期しました.
今回の標本では標本平均がおよそ31,660円となりました.
一方でガンマ分布の平均は形状母数と尺度母数の積で示されるようです.

y = gamma.pdf(x, a=2.5, scale=10000)

このコードに着目すると、それは25,000円であることが分かります.
したがって、今回の標本では、標本平均は母集団の平均から乖離しているということになります. もう少し情報を付け加えると、母集団の平均よりも大きな数値が標本平均として得られました.

ここで再び神の視点になり、標本サイズ23の抽出を何度も行い標本平均の分布を描画してみたいと思います.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib #日本語化matplotlib
import seaborn as sns

sns.set(font="IPAexGothic") #日本語フォント

from scipy.stats import gamma

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8), sharex=True)

x = np.linspace(0, 100000, 1000)
y = gamma.pdf(x, a=2.5, scale=10000)

ax1.set_title('30代の毎月の飲み会代の分布')
ax1.set_ylabel('確率密度関数')
ax2.set_ylabel('標本平均のヒスグラム')

a = 2.5
scale = 10000
gamma_mean = a * scale 
ax1.axvline(gamma_mean, color='r', linestyle='--', label=f'ガンマ分布の平均: {gamma_mean:.2f}')
ax2.axvline(gamma_mean, color='r', linestyle='--', label=f'ガンマ分布の平均: {gamma_mean:.2f}')


# 2.5%点と97.5%点を求める
low = gamma.ppf(0.025, a=2.5, scale=10000)
high = gamma.ppf(0.975, a=2.5, scale=10000)


# 標本をプロットする
np.random.seed(0) # 乱数生成器の状態を固定するために追加した行です。
sample = gamma.rvs(a=2.5, scale=10000, size=23)

means = []
for i in range(1000):
 sample = gamma.rvs(a=2.5, scale=10000, size=23)
 means.append(np.mean(sample))
 
ax2.hist(means)
ax1.plot(x, y)
plt.show()

今回の場合、母集団が従うガンマ分布の最頻値の右側に標本平均の分布が出現するようです.
Bingチャットの上限に達したようなので今回の可視化はここまでです.