面向过程给出《贝叶斯思维：统计建模的Python学习法》——二维彩球问题学习代码

背景

给出读《艾伯特贝叶斯思维：统计建模的Python学习法.pdf》的时候，写的代码，以面向过程的方式给出。
本章彩弹问题，求似然度的时候，假设已知隐藏点时，射手等概率从各个角度射击。

代码

导入常见模块

# %load "E:\桌面space\临时数据\python\个人自定义模块\ImportFile.py"
# Standard Scientific Import
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib.pyplot import plot as plot
import sklearn
import seaborn as sns
import sys
import patsy

# 个人代码测试路径
sys.path.append(r"C:\Users\Administrator\PycharmProjects\QY_TS_Quant")
from QY_plot import *

plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文字体设置-黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
sns.set(font='SimHei',font_scale=1.25,style="ticks",rc={"xtick.major.size": 3, "ytick.major.size": 3})# 解决Seaborn中文显示问题

# Jupyter 默认设置
%matplotlib inline
%config InlineBackend.figure_format="retina"
%config InlineBackend.rc = {"figure.figsize": (7.5,4.5)}

# 多列输出
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

中间函数

from functools import reduce
import operator

# 假设在alpha,beta的射手，随机等概率可以往各个方向设计
def StrafingSpeed(alpha, beta, x):
    return (np.square(beta)+np.square(alpha-x))/beta

def normlize(x):
    return x/np.sum(x)

def MakePmf(alpha, beta, w=30):
    pmf = list(map(lambda v:1./ StrafingSpeed(alpha, beta, v), np.arange(w+1)))
    return normlize(pmf)

def likehood(AX, AY, x, w=30):
    AZS = reduce(operator.add, map(lambda v:1./StrafingSpeed(AX, AY, v),np.arange(w+1)))
    AZ =  1./StrafingSpeed(AX, AY, x)
    return  AZ/AZS

# 围场宽30， 长50
w = 30
h = 50

# 宽边 中弹位置
x = [15, 16, 18, 21]

假定射击点时，墙上各点被击打概率

alpha = 10
beta = 10
plt.plot(MakePmf(alpha, beta))

二维射击点分布

def update(AP,x):
    if not np.isscalar(x):
        for xi in x:
            AP = update(AP, xi)
    else:
        AZ = likehood(AX, AY, x)
        AP = AP*AZ
        AP /= AP.sum()
    return AP

X = np.arange(w+1)
Y = np.arange(1, h+1)
AX, AY = np.meshgrid(X, Y)

AP = np.ones_like(AX)
AP = update(AP, x)

plt.contour(AX,AY,AP,80)
plt.colorbar()

边缘分布

fig, axes = plt.subplots(1, 2, figsize=(18, 6))
axes[0].plot(AX[0], AP.sum(axis=0), label="alpha")
axes[1].plot(AY[:,0], AP.sum(axis=1), label="beta")
axes[0].legend();
axes[1].legend();
plt.suptitle("概率密度");

fig, axes = plt.subplots(1, 2, figsize=(18, 6))
axes[0].plot(AX[0], AP.sum(axis=0).cumsum(), label="alpha")
axes[1].plot(AY[:,0],AP.sum(axis=1).cumsum(), label="beta")
axes[0].legend();
axes[1].legend();
plt.suptitle("分布函数");

def confintv(cdf, p):
    t = []
    flag = 0
    for i,v in enumerate(cdf):
        if flag==0:
            if v>=0.5-p/2:
                t.append(i)
                flag=1
        elif flag==1:
            if v>=0.5+p/2:
                t.append(i)
                break
    return t

# alpha的 50%置信区间
AX[0][confintv(AP.sum(axis=0).cumsum(), 0.5)]

array([14, 21])

# beta的50%置信区间
AY[:,1][confintv(AP.sum(axis=1).cumsum(), 0.5)]

array([ 5, 31])

置信区间

Acolor = np.zeros_like(AP)
s = np.argsort(AP.ravel())
for v in np.searchsorted(AP.flat[s].cumsum(),[0.25, 0.5, 0.75]):
    Acolor.flat[s[:v]] += 0.25

plt.imshow(Acolor, origin=True, aspect=0.5)

clf = plt.contour(AX, AY, Acolor,3);
plt.clabel(clf, fontsize=10);
plt.contourf(AX, AY, Acolor, 3);