統計手法をPythonで実装

二項分布
```Python
# Program Name: binomial_distribution.py
# Creation Date: 20250602
# Overview: Compute and visualize binomial distribution probabilities
# Usage: Run in Python environment with scipy and matplotlib

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

# パラメータ設定 / Set parameters
n = 10         # 試行回数 (number of trials)
p = 0.5        # 成功確率 (probability of success per trial)
k = 3          # 成功回数 (number of successes)

# 確率計算 / Compute P(X = k)
prob = binom.pmf(k, n, p)
print(f"P(X = {k}) when n = {n}, p = {p}:", prob)

# 確率分布の全体表示 / Plot PMF
x = np.arange(0, n+1)
y = binom.pmf(x, n, p)

plt.figure(figsize=(8, 5))
plt.bar(x, y, color='skyblue', edgecolor='black')
plt.title(f'Binomial Distribution: n={n}, p={p}')
plt.xlabel('Number of successes (k)')
plt.ylabel('Probability P(X=k)')
plt.grid(True, axis='y')
plt.show()

```

ポアソン分布
```Python
# Program Name: poisson_defect_probability.py
# Creation Date: 20250602
# Overview: Calculate probability of 1 defect using Poisson distribution
# Usage: Run in Python with scipy installed

from scipy.stats import poisson

# Parameters
n = 10
p = 1 / 200
lam = n * p  # λ = expected number of defects

# Probability that exactly 1 defect is found
prob_1_defect = poisson.pmf(1, lam)
print(f"P(X = 1) = {prob_1_defect:.5f}")

```

幾何分布
```Python
# Program Name: geometric_distribution_dice.py
# Creation Date: 20250602
# Overview: Compute and visualize geometric distribution (dice example)
# Usage: Run in Python with scipy and matplotlib installed

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

# 成功確率（サイコロで1が出る確率）
p = 1/6

# 試行回数の上限（例：1～10回）
k_vals = np.arange(1, 11)

# 幾何分布の確率計算
probs = geom.pmf(k_vals, p)

# 出力確認（例：3回目で初めて成功する確率）
print(f"P(X = 3) = {probs[2]:.5f}")  # インデックス2 → 3回目

# プロット
plt.figure(figsize=(8, 5))
plt.bar(k_vals, probs, color='lightgreen', edgecolor='black')
plt.title('Geometric Distribution: P(X=k) for p=1/6 (Dice)')
plt.xlabel('Trial number k (First success on k-th trial)')
plt.ylabel('Probability P(X=k)')
plt.xticks(k_vals)
plt.grid(axis='y')
plt.show()

```

超幾何分布
```Python
# Program Name: hypergeometric_distribution.py
# Creation Date: 20250602
# Overview: Calculate hypergeometric distribution and apply to examples
# Usage: Run in Python with scipy and matplotlib installed

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

# ✅ 例題1：赤玉と白玉の例
# 全体 N, 赤玉 K, 抽出 n, 成功 k
N = 13      # 8赤 + 5白
K = 8       # 赤玉
n = 3       # 抜き取る数
k = 2       # 赤玉2個が出る確率

prob1 = hypergeom.pmf(k, N, K, n)
print(f"赤玉が2個含まれる確率: {prob1:.4f}")  # ≈ 0.4242

# ✅ 例題2：印付き魚の個体数推定
# 捕獲再捕獲法の逆算（推定）
marked_ratio = 5 / 20   # = 1/4
marked_total = 100
estimated_population = marked_total / marked_ratio
print(f"推定される湖の魚の総数: {estimated_population:.0f}")  # = 400

# ✅ グラフ描画（例題1の分布）
x = np.arange(0, n+1)
y = hypergeom.pmf(x, N, K, n)

plt.figure(figsize=(7, 4))
plt.bar(x, y, color='lightblue', edgecolor='black')
plt.title('Hypergeometric Distribution: P(X=k) for Red Balls')
plt.xlabel('Number of red balls drawn (k)')
plt.ylabel('Probability P(X=k)')
plt.grid(True, axis='y')
plt.show()

```
指数分布

```Python
# 再実行：指数分布の可視化（lambda = 0.1, 0.5, 1, 2, 5, 10）

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

# 時間（hour単位）0～1時間（= 60分）まで
t = np.linspace(0, 1, 500)

# λの例：0.1, 0.5, 1, 2, 5, 10
lambdas = [0.1, 0.5, 1, 2, 5, 10]

plt.figure(figsize=(10, 6))
for lam in lambdas:
    pdf = expon.pdf(t, scale=1/lam)  # scale = 1/λ
    plt.plot(t * 60, pdf, label=f'λ = {lam}')  # x軸は分単位

plt.title('Exponential Distribution (Probability Density)')
plt.xlabel('Time until next event (minutes)')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.show()

```
https://qiita.com/DeepMata/items/45678d3af97c01a3a277

1標本t検定
「平均寿命が80歳という仮説が正しいか？」
調査対象：10人の平均寿命データ
```Python
from scipy.stats import ttest_1samp

# 寿命データ（例）
data = [82, 79, 81, 80, 77, 85, 83, 78, 80, 79]

# t検定（母平均 = 80）
t_stat, p_value = ttest_1samp(data, popmean=80)

print(f"t統計量: {t_stat:.3f}")
print(f"p値: {p_value:.3f}")

```
母分散の検定（カイ二乗分布）
ある製品の重さ（単位：g）にばらつきがあるかを調べたい。
このとき、実際のばらつきが1.5に等しいかを検定せよ。



```Python
import numpy as np
from scipy.stats import chi2

# 標本データ
data = [101.2, 100.8, 101.0, 100.9, 101.1, 100.7, 101.3, 100.6, 101.0, 100.9]
n = len(data)
s2 = np.var(data, ddof=1)          # 標本分散
sigma0_sq = 1.5                    # 仮定する母分散

# カイ二乗統計量
chi2_stat = (n - 1) * s2 / sigma0_sq
df = n - 1                         # 自由度

# 両側検定のp値
p = 2 * min(chi2.cdf(chi2_stat, df), 1 - chi2.cdf(chi2_stat, df))

print(f"標本分散 s² = {s2:.4f}")
print(f"カイ二乗統計量 = {chi2_stat:.4f}")
print(f"自由度 = {df}")
print(f"p値 = {p:.4f}")

```
ある教育法の効果を比較するため、異なる学習方法A・B・Cを3つのクラスに適用し、テストの得点を測定した。
「平均点に有意差があるか」を検定せよ（有意水準5%）

```Python
from scipy.stats import f_oneway

# 各グループのデータ
group_A = [78, 75, 80, 74, 77]
group_B = [85, 90, 78, 82, 89]
group_C = [82, 88, 85, 87, 83]

# ANOVA一元配置
f_stat, p_value = f_oneway(group_A, group_B, group_C)

print(f"F値 = {f_stat:.4f}")
print(f"p値 = {p_value:.4f}")

# 結論
if p_value < 0.05:
    print("⇒ 有意差あり（少なくとも1群に平均差）")
else:
    print("⇒ 有意差なし（平均は統計的に同等）")

```
A/Bテスト：メール開封率、ランディングページ

```python
# Program Name: ab_test_conversion.py
# Creation Date: 20250603
# Overview: Perform A/B test (two-proportion z-test) on email open or conversion rates
# Usage: Execute in Python environment (Jupyter, Colab, or script) with required packages installed

# -------------------- Install Required Libraries --------------------
!pip install numpy scipy matplotlib

# -------------------- Import Libraries --------------------
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# -------------------- Parameter Settings --------------------
# グループA / Group A (Control)
n_A = 1000  # サンプル数 / sample size
x_A = 230   # 成功数（例：開封者数）/ number of successes (e.g., opens)

# グループB / Group B (Treatment)
n_B = 980
x_B = 274

# 有意水準 / Significance level
alpha = 0.05

# -------------------- Compute Proportions --------------------
p_A = x_A / n_A  # グループAの開封率 / proportion of Group A
p_B = x_B / n_B  # グループBの開封率 / proportion of Group B
p_pool = (x_A + x_B) / (n_A + n_B)  # プールされた比率 / pooled proportion

# -------------------- Z検定統計量の計算 / z-test statistic --------------------
# 標準誤差の計算 / standard error
se = np.sqrt(p_pool * (1 - p_pool) * (1/n_A + 1/n_B))

# zスコアの計算 / z-score
z = (p_B - p_A) / se

# 両側検定のp値 / two-sided p-value
p_value = 2 * (1 - norm.cdf(abs(z)))

# -------------------- 結果出力 / Print Result --------------------
print(f"Group A Open Rate (p_A): {p_A:.4f}")
print(f"Group B Open Rate (p_B): {p_B:.4f}")
print(f"Z-statistic: {z:.4f}")
print(f"P-value: {p_value:.4f}")

if p_value < alpha:
    print("結果：有意差あり（BはAと統計的に異なる）")
else:
    print("結果：有意差なし（BはAと統計的に異ならない）")

# -------------------- グラフ描画（検定結果の視覚化） / Plot Test --------------------
x = np.linspace(-4, 4, 1000)
y = norm.pdf(x)

plt.plot(x, y, label="Standard Normal Distribution")
plt.axvline(z, color='red', linestyle='--', label=f"Z = {z:.2f}")
plt.axvline(-z, color='red', linestyle='--')
plt.fill_between(x, y, where=(abs(x) > norm.ppf(1 - alpha/2)), color='gray', alpha=0.3, label='Rejection Region')

plt.title("Z-Test for Two Proportions")
plt.xlabel("Z-value")
plt.ylabel("Density")
plt.legend()
plt.grid(True)
plt.show()

```

売上データからの需要予測（移動平均・指数平滑）
```python
# Program Name: sales_forecast_smoothing.py
# Creation Date: 20250603
# Overview: Demand forecasting using moving average and exponential smoothing
# Usage: Run in Python with pandas, matplotlib installed

# -------------------- Install Required Libraries --------------------
!pip install pandas matplotlib

# -------------------- Import Libraries --------------------
import pandas as pd
import matplotlib.pyplot as plt

# -------------------- Parameter Settings --------------------
# 予測対象データ：仮の売上データ / Sample sales data
data = {
    "Month": pd.date_range(start="2023-01-01", periods=12, freq="M"),
    "Sales": [120, 135, 150, 160, 155, 180, 170, 190, 200, 210, 195, 220]
}
df = pd.DataFrame(data)
df.set_index("Month", inplace=True)

# 移動平均のウィンドウサイズ / Window size for moving average
ma_window = 3

# 指数平滑の平滑化係数 / Smoothing parameter alpha
alpha = 0.3

# -------------------- Moving Average / 移動平均 --------------------
df["Moving_Avg"] = df["Sales"].rolling(window=ma_window).mean()

# -------------------- Exponential Smoothing / 指数平滑法 --------------------
df["Exp_Smooth"] = df["Sales"].ewm(alpha=alpha, adjust=False).mean()

# -------------------- Plot Result / 結果のグラフ表示 --------------------
plt.figure(figsize=(10, 6))
plt.plot(df["Sales"], marker="o", label="Actual Sales")
plt.plot(df["Moving_Avg"], marker="s", linestyle="--", label=f"{ma_window}-Month Moving Average")
plt.plot(df["Exp_Smooth"], marker="^", linestyle="-.", label=f"Exponential Smoothing (α={alpha})")

plt.title("Sales Forecasting with Moving Average and Exponential Smoothing")
plt.xlabel("Month")
plt.ylabel("Sales")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

```

KPI分析とPDCA：統計的モニタリング
```python
# Program Name: kpi_monitoring_pdca.py
# Creation Date: 20250603
# Overview: Statistical KPI monitoring and PDCA cycle evaluation
# Usage: Run in Python with pandas, matplotlib installed

# -------------------- Install Required Libraries --------------------
!pip install pandas matplotlib numpy

# -------------------- Import Libraries --------------------
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# -------------------- Parameter Settings --------------------
# 仮のKPI（コンバージョン率）のデータを生成 / Sample monthly KPI (conversion rate) data
np.random.seed(0)
months = pd.date_range(start="2023-01", periods=12, freq="M")
kpi_values = np.random.normal(loc=0.08, scale=0.01, size=12)  # 平均8%, 標準偏差1%

# 管理限界線（±3σ）/ Control limits (±3σ)
mu = np.mean(kpi_values)
sigma = np.std(kpi_values)
ucl = mu + 3 * sigma  # Upper Control Limit
lcl = mu - 3 * sigma  # Lower Control Limit

# -------------------- Create DataFrame --------------------
df = pd.DataFrame({"Month": months, "ConversionRate": kpi_values})
df.set_index("Month", inplace=True)
df["UCL"] = ucl
df["LCL"] = lcl
df["Mean"] = mu

# -------------------- Plot KPI Chart --------------------
plt.figure(figsize=(10, 6))
plt.plot(df.index, df["ConversionRate"], marker="o", label="Conversion Rate")
plt.plot(df.index, df["UCL"], linestyle="--", color="red", label="UCL (+3σ)")
plt.plot(df.index, df["LCL"], linestyle="--", color="blue", label="LCL (-3σ)")
plt.plot(df.index, df["Mean"], linestyle=":", color="gray", label="Mean")

plt.title("KPI Monitoring Chart: Conversion Rate")
plt.xlabel("Month")
plt.ylabel("Conversion Rate")
plt.ylim(0.04, 0.12)
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# -------------------- PDCA支援コメント / PDCA Action Guidance --------------------
# KPIがUCL/LCLを超えた場合の警告を出力 / Output warning if control limits are breached
alerts = df[(df["ConversionRate"] > ucl) | (df["ConversionRate"] < lcl)]
if not alerts.empty:
    print("⚠️ 以下の月でKPIが管理限界を逸脱しています（改善アクションが必要）:")
    print(alerts[["ConversionRate"]])
else:
    print("✅ 全期間においてKPIは管理限界内にあります（現行施策を維持）")

```