统计_贝叶斯_我的行李箱在哪里(根据《为什么关于因果关系的新科学》中第三章 )

一、$p(我没拿到行李|上飞机)$ 列表

def wait_time(i:int) ->int:
    """
    假定总共等待10分钟
    每等待1分钟，行李在飞机上的概率下降10%
    p(我没拿到行李|上飞机)
    """
    return 1 - i / 10

# p(我没拿到行李|上飞机) list
time_lst = list(range(0, 10))
t_p_lst = [wait_time(i) for i in time_lst]
print(t_p_lst)


"""
[1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.30000000000000004, 0.19999999999999996, 0.09999999999999998]
"""

二、一个时间点的逆概率计算

由代码大师指正:

将 p_npag = t_p * 0.5 + 1 * (1-0.5) 修正为 p_npag = t_p * p_fj + 1 * (1-p_fj)

def pag_not_in_plane(time_i, p_fj = 0.5):
    """
    为什么——关于因果关系的新科学
    一个时间点计算 逆概率
    """
    # 没上飞机概率模型 
    ## p(上飞机|我没拿到行李)p(我没拿到行李) = p(我没拿到行李|上飞机)p(上飞机)
    ## 1. p(我没拿到行李) = p(我没拿到行李|上飞机) + p(我没拿到行李|没上飞机)
    t_p =  wait_time(time_i)
    p_npag = t_p * p_fj + 1 * (1-p_fj) # 修正
    ## 2. p(上飞机) p_fj = 0.5
    ## 3. p(上飞机|我没拿到行李)
    p_taget = p_fj * t_p / p_npag
    return p_taget


p_taget_l = [pag_not_in_plane(i) for i in time_lst]
print([round(i, 2) for i in p_taget_l])

"""
[0.5, 0.47, 0.44, 0.41, 0.37, 0.33, 0.29, 0.23, 0.17, 0.09]
"""

三、绘制等待时期心里预期($p(行李上飞机|我没拿到行李)$)

import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.bar(range(0,10), p_taget_l)
for (x, y) in zip(range(0,10), p_taget_l):
    plt.text(x, y-0.02, round(y,2), ha = 'center')
plt.show()

3.1 再观测心理预期环比变化情况

# 环比变化
import copy
def Hb_value(lst):
    lst1 = [lst[0]]
    lst_t = copy.deepcopy(lst)
    lst_t.append(0)
    lst1.extend(lst_t)
    out_l = [lst_t[i] - lst1[i] for i in range(len(lst_t))]
    out_l.pop()
    return out_l

hb_lst = Hb_value(p_taget_l)
"""
>>> print([round(i, 3) for i in hb_lst])
[0.0, -0.026, -0.029, -0.033, -0.037, -0.042, -0.048, -0.055, -0.064, -0.076]
"""


def plot_rate_div(time_lst, p_taget_l): 
    # 环比值
    hb_lst = Hb_value(p_taget_l)
    # 画图
    fig, ax1 = plt.subplots(figsize = (16, 8))

    # 左轴
    ax1.bar(time_lst, p_taget_l, color='steelblue', alpha=0.5, label = 'rate')
    for (x, y) in zip(time_lst, p_taget_l):
        ax1.text(x, y-0.02, round(y,2), ha = 'center')
    ax1.set_xlabel('time')
    ax1.set_ylabel('rate')
    plt.legend(loc='upper left')

    # 右轴
    ax2 = ax1.twinx()
    ax2.plot(time_lst[1:], hb_lst[1:], '-or', label = 'div_values')
    ax2.set_ylabel('div rate')
    # ax2.set_ylim(-0.2, 0.1)
    plt.legend(loc='upper right')
    plt.show()

可以发现心里预期前期下降慢后期下降快，和我们的一般心理预期相符合

四、如果更改行李箱上飞机的概率(0.8)的话

p_taget_l_cg = [pag_not_in_plane(i, 0.8) for i in time_lst]
plot_rate_div(time_lst, p_taget_l_cg)

由此可见我们可以等7分钟左右，如果7分钟还没有出现，那么心理预期会大幅下降。