贝叶斯统计 最大后验密度可信集 91,95页
import numpy as np
import sympy as sy
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
alpha_bar=1-0.05
x=np.arange(80,140,0.01)
x=np.around(x,decimals=2)
#pdf
y=np.around(1/(np.sqrt(2*np.pi)*8.32)*np.exp(-1*((x-110.38)**2)/(2*8.32**2)),decimals=9)#N(110.38,8.32^2)
plt.plot(x,y)
#符号变量
sy_x=sy.symbols('x')
sy_y=1/(sy.sqrt(2*sy.pi)*8.32)*sy.exp(-1*((sy_x-110.38)**2)/(2*8.32**2))#N(110.38,8.32^2)
#初始赋值
k=np.median(y)
theta_1_k=x[np.argwhere(y==k)[0]]
theta_2_k=x[np.argwhere(y==k)[1]]
s=sy.integrate(sy_y,(sy_x,theta_1_k,theta_2_k))
n=1
while True:
if s-alpha_bar>0.0005:
k=min(y[y>k])
print(n,":",s-alpha_bar)
n=n+1
elif alpha_bar-s>0.0005:
k=max(y[y
结果：
迭代128次 可信区间[94.1，126.66]
