java 后验条件_基于贝叶斯后验优化的曲线拟合

我试图使用Python PYMC3包在我的数据上创建后验预测分布，得到累积和条件概率作为最终结果 .

我正在研究3种预期寿命：12个月，24个月和36个月 . 个人预期寿命群体在其生命中发生死亡时具有不同的历史形态 .

例如，这是基于历史信息的24个月预期寿命模式：

所以我正在探索用于曲线拟合的贝叶斯方法，并且一直在尝试使用负二项分布来创建适合这些数据的曲线 . (我认为lognormal会更合适，但我没有机会调整我的代码 .

下面是我用来拟合曲线的代码/逻辑：

# life_expectancy = 12, 24, 36

# dead = 1, 0

indiv_traces = {}

# Convert categorical variables to integer

le = preprocessing.LabelEncoder()

participants_idx = le.fit_transform(df_comb_clean[(df_comb_clean['dead']==1)]['life_expectancy'])

participants = le.classes_

n_participants = len(participants)

for p in participants:

with pm.Model() as model:

alpha = pm.Uniform('alpha', lower=0, upper=100)

mu = pm.Uniform('mu', lower=0, upper=100)

data = df_comb_clean[(df_comb_clean['dead']==1) & (df_comb_clean['life_expectancy']==p)]['month'].values

y_est = pm.NegativeBinomial('y_est', mu=mu, alpha=alpha, observed=data)

y_pred = pm.NegativeBinomial('y_pred', mu=mu, alpha=alpha)

start = pm.find_MAP()

step = pm.Metropolis()

trace = pm.sample(20000, step, start=start, progressbar=True)

indiv_traces[p] = trace

结果：

# Optimization life_expectancyinated successfully.

# Current function value: 133.738083

# Iterations: 11

# Function evaluations: 19

# Gradient evaluations: 19

# 100%|██████████| 20000/20000 [00:12<00:00, 1659.15it/s]

# Optimization life_expectancyinated successfully.

# Current function value: 679.448463

# Iterations: 17

# Function evaluations: 25

# Gradient evaluations: 25

# 100%|██████████| 20000/20000 [00:12<00:00, 1614.50it/s]

# Optimization life_expectancyinated successfully.

# Current function value: 812.109878

# Iterations: 18

# Function evaluations: 201

# Gradient evaluations: 198

# 100%|██████████| 20000/20000 [00:12<00:00, 1570.42it/s]

现在我正在绘制我的后验预测分布图：

combined_y_pred = np.concatenate([v.get_values('y_pred') for k, v in indiv_traces.items()])

x_lim = 24

y_pred = trace.get_values('y_pred')

fig = plt.figure(figsize=(12,6))

fig.add_subplot(211)

fig.add_subplot(211)

_ = plt.hist(combined_y_pred, range=[5, x_lim], bins=x_lim, histtype='stepfilled', color=colors[1])

_ = plt.xlim(4, x_lim)

_ = plt.ylim(0, 5000)

_ = plt.ylabel('Frequency')

_ = plt.title('Posterior predictive distribution')

fig.add_subplot(212)

# ter

# df_comb_co['month'].values,

_ = plt.hist(df_comb_clean[df_comb_clean['dead']==1]['month'].values,range=[5, x_lim], bins=x_lim, histtype='stepfilled')

_ = plt.xlim(4, x_lim)

_ = plt.xlabel('Month')

_ = plt.ylim(0, 50)

_ = plt.ylabel('Frequency')

_ = plt.title('Distribution of observed data')

plt.tight_layout()

然后得到以下输出：

所以现在，我想提取我的结果并将它们转换为适合我的初始数据的条件(基于月)曲线 . 我通过以下代码进行了基本尝试，预期寿命为24个月：

def life_expectancy_y_pred(life_expectancy):

"""Return posterior predictive for person"""

ix = np.where(participants == life_expectancy)[0][0]

return trace['y_pred']

life_expectancy = 24

x = np.linspace(4, life_expectancy, num=life_expectancy)

num_samples = float(len(life_expectancy_y_pred(life_expectancy)))

prob_lt_cum_x = [sum(life_expectancy_y_pred(life_expectancy) < i)/num_samples for i in x]

下面是我的结果，其中蓝色是实际的，橙色是负二项式拟合我在python中手动完成，黄色是我从贝叶斯优化过程得到的 .

请告诉我我做错了什么，因为我对贝叶斯的适应性非常糟糕 . 我希望右尾相对于它的位置，但是尖端类似于橙色线 .

我不想放弃这个过程 .