独立事件组合概率与收益估算：构建总收益概率分布曲线

jobs = ['job1', 'job2', 'job3', 'job4', 'job5']
probabilities = [0.1, 0.1, 0.4, 0.6, 0.2] # 对应每个项目的成功概率
hours = [1, 10, 43, 2, 5] # 对应每个项目成功后的工时收益
min_hours_desired = 10 # 我们希望计算达到或超过这个工时的概率

scenarios = []
jobs_len = len(jobs)
for i in range(2**jobs_len):
    # 将整数i转换为二进制字符串，并用0填充到jobs_len长度
    scenario = bin(i).split('b')[1].zfill(jobs_len)
    scenarios.append(scenario)

print(f"生成的总情景数: {len(scenarios)}")
# 示例：打印前几个情景
# for s in scenarios[:5]:
#     print(s)

scenario_outcomes = []
for scenario in scenarios:
    scenario_hours_won = 0
    scenario_probability = 1.0 # 使用浮点数确保精确计算
    for j, b in enumerate(scenario):
        if b == '0': # 项目失败
            scenario_probability *= (1 - probabilities[j])
        else: # 项目成功
            scenario_probability *= probabilities[j]
            scenario_hours_won += hours[j]
    scenario_outcomes.append((scenario, scenario_probability, scenario_hours_won))

# 打印部分情景的计算结果作为示例
print("\n部分情景的概率与收益：")
for outcome in scenario_outcomes[:5]:
    print(f"情景: {outcome[0]}, 概率: {outcome[1]:.6f}, 总工时: {outcome[2]}")

prob_desired_hours = sum([o[1] for o in scenario_outcomes if o[2] > min_hours_desired])
print(f'\n获得大于 {min_hours_desired} 小时的总概率: {prob_desired_hours:.6f}')

# 验证所有情景的概率之和是否为1
prob_check = sum([o[1] for o in scenario_outcomes])
print(f'所有情景概率之和 (应为1): {prob_check:.6f}')

import collections

# 使用defaultdict来累加相同工时的概率
payout_probabilities = collections.defaultdict(float)

for _, prob, payout in scenario_outcomes:
    payout_probabilities[payout] += prob

# 将结果转换为字典并按工时排序，方便后续绘图
sorted_payout_probabilities = dict(sorted(payout_probabilities.items()))

import json
print("\n总工时及其对应概率分布：")
print(json.dumps(sorted_payout_probabilities, indent=2))

# 进一步，我们可以用matplotlib等库将这个分布绘制成柱状图或散点图。
# import matplotlib.pyplot as plt
#
# payouts = list(sorted_payout_probabilities.keys())
# probabilities_values = list(sorted_payout_probabilities.values())
#
# plt.figure(figsize=(10, 6))
# plt.bar(payouts, probabilities_values, width=1.0, edgecolor='black')
# plt.xlabel('总工时')
# plt.ylabel('发生概率')
# plt.title('总工时概率分布')
# plt.xticks(rotation=45)
# plt.grid(axis='y', linestyle='--', alpha=0.7)
# plt.tight_layout()
# plt.show()