Think Bayes Note 12: 彩弹

Author: nex3z 2018-08-28

1. 问题描述

　　在彩弹射击游戏中，参加游戏的双方使用装有彩弹的玩具枪相互射击，彩弹在命中目标后，会在目标上染色。假设游戏场景是一个 30 英尺长、50 英尺宽的房间，你发现 30 英尺长一面的墙上，在距离墙角 15 英尺、16 英尺、18 英尺和 21 英尺的位置上各有一个彩弹的痕迹，据此推测对手的位置。

2. 场景分析

　　问题场景如下图所示。

假设墙角为原点，射击方的位置为 $(\alpha, \beta)$，射击点为 $x$，射击角度为 $\theta$。由图可知：

可得设计角度 $\theta$ 为：

　　另外，式 (1) 等号两边对 $\theta$ 求导，可得：

式 (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()

图像为：

可见对于更大的置信区间，分布更偏向右侧。