切断正規分布のパラメータ計算

与えられた期待値と分散から切断正規分布のパラメータを求める方法です。

scipyのoptimizeを使います。

準備1 公式に従って関数を用意

import math
from scipy import stats, optimize
def dist_E(lower, upper, mu, sigma):
  x0 = (lower - mu) / sigma
  x1 = (upper - mu) / sigma
  return mu + (stats.norm.pdf(x0) - stats.norm.pdf(x1)) / (stats.norm.cdf(x1) - stats.norm.cdf(x0)) * sigma

def dist_V(lower, upper, mu, sigma):
  x0 = (lower - mu) / sigma
  x1 = (upper - mu) / sigma
  term1 = (x0 * stats.norm.pdf(x0) - x1 * stats.norm.pdf(x1)) / (stats.norm.cdf(x1) - stats.norm.cdf(x0))
  term2 =  (stats.norm.pdf(x0) - stats.norm.pdf(x1)) / (stats.norm.cdf(x1) - stats.norm.cdf(x0))
  return sigma * sigma * (1.0 + term1 - term2 * term2)

準備2
推定条件を準備

a= 1 分布の下限
b= 5 分布の上限
E= 4.00 分布の期待値
V= 0.5 分布の分散

準備3
目的関数を準備

def fun(x):
  return [dist_E(a, b, x[0], x[1]) - E, dist_V(a, b, x[0], x[1]) -V]

パラメータ計算

sol = optimize.root(fun, x0=[E, math.sqrt(V)])
x = sol.x
mu = x[0]
sigma = x[1]
print(x)
print(dist_E(a, b, mu, sigma))
print(dist_V(a, b, mu, sigma))

結果

[4.61408088 1.05919093]
4.0
0.4999999999999998

グラフ

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)

fix_a = (a - mu) / sigma
fix_b = (b - mu) / sigma
print(stats.truncnorm.mean(fix_a, fix_b, loc=mu, scale=sigma,  moments='mvsk'))
print(stats.truncnorm.var(fix_a, fix_b, loc=mu, scale=sigma,  moments='mvsk'))
x = np.linspace(stats.truncnorm.ppf(0.01, fix_a, fix_b, loc=mu, scale=sigma), stats.truncnorm.ppf(0.99, fix_a, fix_b, loc=mu, scale=sigma), 100)
ax.plot(x, stats.truncnorm.pdf(x, fix_a, fix_b, loc=mu, scale=sigma), 'r-', lw=5, alpha=0.6, label='truncnorm pdf')
plt.show()

グラフ結果

4.0
0.4999999999999998