\begin{align*}
p(\mu|D)&\propto  p(\mu)p(D|\mu)\\
&=p(\mu)\prod_{x_i\in D} p(x_i|\mu)\\
&\propto 
\exp\left\{-\frac{(\mu-\mu_0)^2}{2\sigma_0^2}
-\sum_i\frac{(x_i-\mu)^2}{2\sigma^2}\right\}\\
&\propto \exp\left\{-\frac{1}{2}\left((n/\sigma^2+1/\sigma_0^2)\mu^2-2(\sum_i x_i/\sigma^2+\mu_0/\sigma_0^2)\mu+\mathrm{const.}\right)\right\}\\
&\propto \exp\left\{-\frac{1}{2}(n/\sigma^2+1/\sigma_0^2)\left(\mu-\frac{\sum_i x_i/\sigma^2+\mu_0/\sigma_0^2}{n/\sigma^2+1/\sigma_0^2}\right)^2\right\}
\end{align*}

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, ttest_ind

# 実験設定
num_a = 200
num_b = 190

true_mu_a = 2
true_mu_b = 3
sigma = 4 # 標準偏差(既知)

data_a=np.random.normal(true_mu_a, sigma, size=num_a)
data_b=np.random.normal(true_mu_b, sigma, size=num_b)

mean_a = data_a.mean()
mean_b = data_b.mean()
print(f"{mean_a=}")
print(f"{mean_b=}")

# ABテスト
test=ttest_ind(data_a, data_b, equal_var=False, alternative="two-sided")
print(f"{test.pvalue=}(two-sided)")

test=ttest_ind(data_a, data_b, equal_var=False, alternative="less")
print(f"{test.pvalue=}(less)")

# ベイジアンAB
sigma0 = sigma
mu0 = 2.5

x = np.linspace(0, 5, 1000)

mu_a = (num_a*mean_a/sigma**2 + mu0/sigma0**2)/(num_a/sigma**2 + 1/sigma0**2)
var_a = 1/(num_a/sigma0**2 + 1/sigma0**2)
sigma_a = np.sqrt(var_a)
mu_b = (num_b*mean_b/sigma**2 + mu0/sigma0**2)/(num_b/sigma**2 + 1/sigma0**2)
var_b = 1/(num_b/sigma0**2 + 1/sigma0**2)
sigma_b = np.sqrt(var_b)

y_a = norm.pdf(x, mu_a, sigma_a)
y_b = norm.pdf(x, mu_b, sigma_b)

plt.plot(x, y_a, label="A")
plt.plot(x, y_b, "--", label="B")
plt.legend()
plt.savefig("normal.png")

# 勝率計算
mu_ab=mu_a-mu_b
sigma_ab=np.sqrt(var_a+var_b)

a_less_than_b_prob=norm.cdf(0, mu_ab, sigma_ab)
print(f"p(A<B) = {a_less_than_b_prob}")

# サンプリングによる勝率計算
simulation_size=100000
x_a = np.random.normal(mu_a, sigma_a, size=simulation_size)
x_b = np.random.normal(mu_b, sigma_b, size=simulation_size)
print(f"p(A<B) = {np.mean(x_a<x_b)}")

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.proportion import proportions_ztest
from scipy.stats import beta

num_a = 200
num_b = 190

true_mu_a = 0.2
true_mu_b = 0.3

data_a=np.random.binomial(n=1, p=true_mu_a, size=num_a)
data_b=np.random.binomial(n=1, p=true_mu_b, size=num_b)

mean_a = data_a.mean()
mean_b = data_b.mean()
print(f"{mean_a=}")
print(f"{mean_b=}")

sum_a = data_a.sum()
sum_b = data_b.sum()

# ABテスト
z_score, p_value = proportions_ztest(
    [sum_a, sum_b], [num_a, num_b], 
    alternative="two-sided"
)
print(f"{p_value=}(two-sided)")

z_score, p_value = proportions_ztest(
    [sum_a, sum_b], [num_a, num_b], 
    alternative="smaller"
)
print(f"{p_value=}(smaller)")

# ベイジアンAB
x = np.linspace(0, 0.5, 1000)
y_a = beta.pdf(x, 1+sum_a, 1+num_a-sum_a)
y_b = beta.pdf(x, 1+sum_b, 1+num_b-sum_b)

plt.plot(x, y_a, label="A")
plt.plot(x, y_b, "--", label="B")
plt.legend()
plt.savefig("beta.png")

# サンプリングによる勝率計算
simulation_size=100000
x_a = np.random.beta(1+sum_a, 1+num_a-sum_a, size=simulation_size)
x_b = np.random.beta(1+sum_b, 1+num_b-sum_b, size=simulation_size)
print(f"p(A<B) = {np.mean(x_a<x_b)}")