Think Bayes Note 12: 彩弹
1. 问题描述
在彩弹射击游戏中，参加游戏的双方使用装有彩弹的玩具枪相互射击，彩弹在命中目标后，会在目标上染色。假设游戏场景是一个 30 英尺长、50 英尺宽的房间，你发现 30 英尺长一面的墙上，在距离墙角 15 英尺、16 英尺、18 英尺和 21 英尺的位置上各有一个彩弹的痕迹，据此推测对手的位置。
2. 场景分析
问题场景如下图所示。
假设墙角为原点，射击方的位置为 $(\alpha, \beta)$，射击点为 $x$，射击角度为 $\theta$。由图可知：
\begin{equation}
\frac{x – \alpha}{\beta} = \tan(\theta) \tag{1}
\end{equation}
可得设计角度 $\theta$ 为：
\begin{equation}
\theta = \tan^{-1}(\frac{x – \alpha}{\beta}) \tag{2}
\end{equation}
另外，式 (1) 等号两边对 $\theta$ 求导，可得：
\begin{equation}
\frac{dx}{d\theta} = \frac{\beta}{\cos^2\theta} \tag{3}
\end{equation}
式 (3) 中的 $\frac{dx}{d\theta}$ 是射击位置 $x$ 随射击角度 $\theta$ 变化的速度，可以称之为“扫射速度”，命中墙上一个特定点的概率与扫射速度负相关。
3. 先验概率
认为射中墙上任意位置的先验概率都相同，定义 PaintBall 如下：
class PaintBall(Suite, Joint):
def __init__(self, alphas, betas, locations):
self.locations = locations
pairs = [(alpha, beta) for alpha in alphas for beta in betas]
Suite.__init__(self, pairs)
...
因为射击位置由 $\alpha$ 和 $\beta$ 两个随机变量共同决定，这里 PaintBall 还继承了表示联合分布的 Joint。
4. 射中任意位置的概率
通过式 (2) 可以由射击位置 $\alpha$、 $\beta$ 和命中位置 $x$ 计算得到射击角度 $\theta$，通过式 (3) 可以进一步得到扫射速度：
def strafing_speed(alpha, beta, x):
theta = math.atan2(x - alpha, beta)
speed = beta / math.cos(theta) ** 2
return speed
由于命中墙上一个特定点的概率与扫射速度负相关，可以据此得到命中任意位置的概率的 PMF：
def make_location_pmf(alpha, beta, locations, name=''):
probs = [1.0 / strafing_speed(alpha, beta, x) for x in locations]
pmf = Pmf(values=locations, probs=probs, name=name)
return pmf
绘制 $\alpha = 10$，$\beta = 10,20,40$ 时命中任意点数的 PMF 如下：
locations = range(31) pmf_beta_10 = make_location_pmf(10, 10, locations, "beta = 10") pmf_beta_20 = make_location_pmf(10, 20, locations, "beta = 20") pmf_beta_40 = make_location_pmf(10, 40, locations, "beta = 40") pmf_beta_10.plot_with([pmf_beta_20, pmf_beta_40])
图像为：
可见最可能命中的位置为 $x = 10$，随着 $\beta$ 的增大，PMF 的范围也逐渐扩大。
5. 似然度
得到命中任意位置的 PMF 后，可以定义似然度函数如下：
class PaintBall(Suite, Joint):
...
def likelihood(self, data, hypo):
alpha, beta = hypo
x = data
pmf = make_location_pmf(alpha, beta, self.locations)
like = pmf.prob(x)
return like
其中 hypo 是射击方位置的假设，data 为命中位置，通过 make_location_pmf() 得到命中任意位置的 PMF 后，求得命中指定位置的概率。
使用以下代码模拟问题中描述的过程：
alphas = range(31) betas = range(1, 51) locations = range(0, 31) paint_ball = PaintBall(alphas, betas, locations) paint_ball.update_set([15, 16, 18, 21])
规定了射击位置 alpha 和 beta 的范围，以及命中位置 locations 的范围后，构造 PaintBall 实例，通过 update_set() 更新命中位置，计算后验概率。
6. 边缘分布
然后可以计算 $\alpha$ 和 $\beta$ 的边缘概率分布 CDF，计算置信区间并绘制图像：
pmf_alpha = paint_ball.marginal(0)
pmf_beta = paint_ball.marginal(1)
print("alpha CI: {}, beta CI: {}"
.format(pmf_alpha.credible_interval(50), pmf_beta.credible_interval(50)))
cdf_alpha = pmf_alpha.make_cdf(name="alpha")
cdf_beta = pmf_beta.make_cdf(name="beta")
cdf_alpha.plot_with([cdf_beta])
其中 marginal() 的定义为：
class Joint(Pmf):
def marginal(self, i, name=''):
pmf = Pmf(name=name)
for vs, prob in self.iter_items():
pmf.incr(vs[i], prob)
return pmf
...
即遍历联合分布中的第 i 个变量在另一个变量取任意值时的概率，构造 PMF。
输出为：
alpha CI: (14.0, 21.0), beta CI: (5.0, 31.0)
图像为：
可见 $\alpha$ 主要集中在 18 附近；$\beta$ 更多地在 10 以内，在 10 以外基本上均匀分布。
7. 条件分布
绘制 \beta = 10, 20, 40 时 \alpha 的条件分布 PMF 如下：
pmf_cond_beta_10 = paint_ball.conditional(0, 1, 10, name="beta = 10") pmf_cond_beta_20 = paint_ball.conditional(0, 1, 20, name="beta = 20") pmf_cond_beta_40 = paint_ball.conditional(0, 1, 40, name="beta = 40") pmf_cond_beta_10.plot_with([pmf_cond_beta_20, pmf_cond_beta_40])
其中 conditional() 的定义为：
class Joint(Pmf):
...
def conditional(self, i, j, val, name=''):
pmf = Pmf(name=name)
for vs, prob in self.iter_items():
if vs[j] != val:
continue
pmf.incr(vs[i], prob)
pmf.normalize()
return pmf
即遍历连个分布中第 i 个变量在第 j 个变量取值为 val 时的概率，构造 PMF。
图像为：
可见 $\beta$ 越小，$\alpha$ 的分布越集中。
8. 置信区间
绘制 25%、50%、75% 置信区间如下：
x = np.arange(0, 31, 1)
y = np.arange(0, 51, 1)
X, Y = np.meshgrid(x, y)
Z = np.zeros(X.shape)
for x, y in interval_25:
Z[y, x] += 1
for x, y in interval_50:
Z[y, x] += 1
for x, y in interval_75:
Z[y, x] += 1
fig = plt.contourf(X, Y, Z, 3, cmap=plt.cm.Blues, antialiased=False)
cbar = plt.colorbar(fig)
cbar.ax.set_yticklabels(['1', '0.75', '0.5', '0.25'])
plt.show()
图像为：
可见对于更大的置信区间，分布更偏向右侧。