《贝叶斯统计》最大后验密度可信集(HPD)Python实现

贝叶斯统计 最大后验密度可信集 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<k])
        print(n,":",alpha_bar-s)
        n=n+1
    else:
        break
    print(theta_1_k,theta_2_k)
    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))

print(theta_1_k,theta_2_k)

结果：

迭代128次 可信区间[94.1，126.66]