負の二項分布のMとK インタラクティブノートブック

https://gyazo.com/91654aaf38056945e4450e9e31d7c5e7

code:python

import numpy as np

import matplotlib.pyplot as plt

from ipywidgets import interact, FloatSlider, IntSlider

def sample_nb(mean, k, size):

"""Sample from a Negative Binomial with given mean and shape k."""

p = k / (k + mean) # success probability

r = k # shape (number of failures)

return np.random.negative_binomial(r, p, size=size)

@interact(

mean_A=FloatSlider(value=5.0, min=0.5, max=50.0, step=0.5, description='Mean A'),

mean_B=FloatSlider(value=5.0, min=0.5, max=50.0, step=0.5, description='Mean B'),

K_A=FloatSlider(value=1.0, min=0.5, max=50.0, step=0.5, description='K A'),

K_B=FloatSlider(value=5.0, min=0.5, max=50.0, step=0.5, description='K B'),

threshold=IntSlider(value=15, min=1, max=150, step=1, description='Threshold'),

)

def play(mean_A, mean_B, K_A, K_B, threshold):

np.random.seed(42)

n_samples = 10_000 # per period per company

periods = 3

samples_A = sample_nb(mean_A, K_A, size=(periods, n_samples))

samples_B = sample_nb(mean_B, K_B, size=(periods, n_samples))

total_A = samples_A.sum(axis=0)

total_B = samples_B.sum(axis=0)

# Visualise distributions

fig, ax = plt.subplots()

bins = np.arange(0, max(total_A.max(), total_B.max()) + 1)

ax.hist(total_A, bins=bins, alpha=0.6, label=f'Company A (M={mean_A}, K={K_A})', density=True)

ax.hist(total_B, bins=bins, alpha=0.6, label=f'Company B (M={mean_B}, K={K_B})', density=True)

ax.axvline(threshold, linestyle='--', label='Threshold')

ax.set_xlabel('Total usage events over 3 periods')

ax.set_ylabel('Density')

ax.set_title('Distribution of summed usage – 3 periods')

ax.legend()

plt.show()

def stats(x):

return {

'mean': np.mean(x),

'std': np.std(x),

'95th pct': np.percentile(x, 95),

f'P(>{threshold})': np.mean(x > threshold),

}

print('Company A stats:', stats(total_A))

print('Company B stats:', stats(total_B))