运筹学练习Python精解——决策论

练习1

某公司有10万元余外资金。如用于开发某个项目估量成功率为95%，成功时一年可获利15%，但一旦失败，有全部丧失资金的危险。如把资金存放到银行中，则可稳得年利4%。为获得更多的信息，该公司求助于咨询公司，咨询费为800元，但咨询意见只是提供参考。拒过去咨询公司类似200例咨询意见实施结果如下表所示，试用决策树法分析：(1)该公司是否值得求助咨询公司；(2)该公司余外资金该如何使用？

咨询意见	投资成功	投资失败	合计
能够投资	150 次	6 次	156 次
不宜投资	22 次	22 次	44 次
合计	172	28	200 次

# 初始概率
P_S = 0.95  # 项目成功概率
P_F = 0.05  # 项目失败概率

# 咨询公司给出的意见概率
P_A_given_S = 150 / 172
P_A_given_F = 6 / 28
P_N_given_S = 22 / 172
P_N_given_F = 22 / 28

# 打印咨询公司给出的意见概率矩阵
print("咨询公司给出的意见概率矩阵:")
print("P(A|S): {:.2f}".format(P_A_given_S))
print("P(A|F): {:.2f}".format(P_A_given_F))
print("P(N|S): {:.2f}".format(P_N_given_S))
print("P(N|F): {:.2f}".format(P_N_given_F))

# 咨询意见的全概率
P_A = P_A_given_S * P_S + P_A_given_F * P_F
P_N = P_N_given_S * P_S + P_N_given_F * P_F

print("P_A: {:.2f}".format(P_A))
print("P_N: {:.2f}".format(P_N))  
      
P_S_given_A = (P_A_given_S * P_S) / P_A
P_F_given_A = (P_A_given_F * P_F) / P_A
P_S_given_N = (P_N_given_S * P_S) / P_N
P_F_given_N = (P_N_given_F * P_F) / P_N

# 打印后验概率矩阵
print("后验概率矩阵:")
print("P(S|A): {:.2f}".format(P_S_given_A))
print("P(F|A): {:.2f}".format(P_F_given_A))
print("P(S|N): {:.2f}".format(P_S_given_N))
print("P(F|N): {:.2f}".format(P_F_given_N))

# 定义收益
gain_S = 15  # 项目成功，收益
loss_F = -100  # 项目失败，损失
gain_bank = 4  # 银行存款收益
cost_of_report = 0.8  # 咨询费用

# 不进行咨询的期望收益
expected_no_report = P_S * gain_S + P_F * loss_F

# 进行咨询后的期望收益
expected_with_report_invest = P_S_given_A * gain_S + P_F_given_A * loss_F 
expected_with_report_not_invest = gain_bank

expected_with_report = P_A * max(expected_with_report_invest, expected_with_report_not_invest) + \
                       P_N * gain_bank- cost_of_report

print("不进行咨询的期望收益: {:.2f}".format(expected_no_report))
print("进行咨询后的期望收益: {:.2f}".format(expected_with_report))

# 比较期望收益，选择最优决策
if expected_with_report > expected_no_report:
    print("最佳决策是: 进行咨询")
else:
    print("最佳决策是: 不进行咨询")

#https://wenku.baidu.com/view/32c8f4f5590216fc700abb68a98271fe900eaf48.html?_wkts_=1718502540796&bdQuery=%E5%86%B3%E7%AD%96%E8%AE%BA+%E9%A3%8E%E9%99%A9%E5%9E%8B%E5%86%B3%E7%AD%96+%E4%B8%8D%E5%AE%8C%E5%85%A8%E4%BF%A1%E6%81%AF+%E4%B9%A0%E9%A2%98+%E4%BE%8B%E9%A2%98
咨询公司给出的意见概率矩阵:
P(A|S): 0.87
P(A|F): 0.21
P(N|S): 0.13
P(N|F): 0.79
P_A: 0.84
P_N: 0.16
后验概率矩阵:
P(S|A): 0.99
P(F|A): 0.01
P(S|N): 0.76
P(F|N): 0.24
不进行咨询的期望收益: 9.25
进行咨询后的期望收益: 11.20
最佳决策是: 进行咨询

练习2

某工程项目按合同应在三个月内完工，其施工费用与工程完工期有关。假定天气是影响能否按期完工的决定因素，如果天气好，工程能按时完工，获利5万元；如果天气不好，不能按时完工，施工单位将被罚款1万元；若不施工就要付出窝工费2千元。根据过去的经验，在计划实施工期天气好的可能性为30%。为了更好地掌握天气情况，可以申请气象中心进行天气预报，并提供同一时期天气预报资料，但需要支付资料费800元。从提供的资料中可知，气象中心对好天气预报准确性为80%，对坏天气预报准确性为90%。问如何进行决策。

解题过程：

先验分析
根据已有资料做出决策损益表。

	施工（d1）	不施工(d2)
好天气 $\theta_1$ (0.3)	5	-0.2
坏天气 $\theta_2$ (0.7)	-1	-0.2
$\mathbf{E}(\mathbf{dj})$	0.8	-0.2

根据期望值准则选择施工方案有利, 相应最大期望收益值 $\mathbf{EMV}$ *(先) $=0.8$

后验概率
气象中心将提供预报此时期内两种天气状态 $x_1$ (好天气)、 $x_2$ (坏天气)将会出现哪一种状态。

先验概率		条件概率		$P(\mathrm{xi} \cap \theta \mathrm{j})$		后验概率
	$P(\theta \mathrm{j})$	$\mathrm{x}1$	$\mathrm{x}2$	$\mathrm{x}1$	$\mathrm{x}2$	$\mathrm{x}1$	$\mathrm{x}2$
好天气 $\theta 1$	0.3	0.8	0.2	0.24	0.06	0.77	0.09
坏天气 $\theta 2$	0.7	0.1	0.9	0.07	0.63	0.23	0.91
				$P(\mathrm{x}1) = 0.31$	$P(\mathrm{x} 2) = 0.69$

后验决策
若气象中心预报天气好（ $x_1$ ）, 则每个方案的最大期望收益值

\begin{aligned} E (d 1 / x 1) = 0.77 \times 5 + 0.23 \times (- 1) = 3.62 \\ E (d 2 / x 1) = 0.77 \times (- 0.2) + 0.23 \times (- 0.2) = - 0.2 \end{aligned}

$\begin{aligned} & E(d1 / x 1)=0.77 \times 5+0.23 \times(-1)=3.62 \\ & E(d2 / x 1)=0.77 \times(-0.2)+0.23 \times(-0.2)=-0.2 \end{aligned}$

若气象中心预报天气不好（ $x_2$ ），各方案的最大期望收益值

\begin{aligned} E (d 1 / x_{2}) = 0.09 \times 5 + 0.91 \times (- 1) = - 0.46 \\ E (d 2 / x_{2}) = 0.09 \times (- 0.2) + 0.91 \times (- 0.2) = - 0.2 \end{aligned}

$\begin{aligned} & E(d1/x_2)=0.09\times 5+0.91 \times(-1)=-0.46 \\ & E(d2/x_2)=0.09 \times(-0.2)+0.91 \times(-0.2)=-0.2 \end{aligned}$

选择 $d1$ 即施工的方案，相应在预报 $x_1$ 时的最大期望收益值 $E(x_1)）=3.62$
选择 $d2$ 即不施工的方案，相应在预报 $x_2$ 时的最大期望收益值 $E(x_2）=-0.2$

计算补充信息的价值：

得到天气预报的情况下，后验决策的最大期望收益值：

EMV*(后) = P (x_{i}) \cdot E (x_{i}) + P (x_{2}) \cdot E (x_{2}) = 0.31 \times 3.62 + 0.69 \times (- 0.2) = 0.9842

$\text{EMV*(后)} = P(x_i) \cdot E(x_i) + P(x_2) \cdot E(x_2) = 0.31 \times 3.62 + 0.69 \times (-0.2) = 0.9842$

则补充的信息价值为：

EMV*(后) - EMV*(先) = 0.9842 - 0.08 = 0.9042

$\text{EMV*(后)} - \text{EMV*(先)} = 0.9842 - 0.08 = 0.9042$

补充信息价值大于信息费（800元），即这种费用是合算的。

import numpy as np

# 定义初始概率
P = np.array([0.3, 0.7])  # P_G, P_B

# 定义气象中心的预测准确性
P_A_given = np.array([0.8, 0.1])  # P_A_given_G, P_A_given_B
P_N_given = np.array([0.2, 0.9])  # P_N_given_G, P_N_given_B

# 计算全概率
P_A = np.dot(P_A_given, P)
P_N = np.dot(P_N_given, P)

# 使用贝叶斯定理计算后验概率
P_G_given_A = (P_A_given[0] * P[0]) / P_A
P_B_given_A = (P_A_given[1] * P[1]) / P_A

P_G_given_N = (P_N_given[0] * P[0]) / P_N
P_B_given_N = (P_N_given[1] * P[1]) / P_N

# 定义收益
gain_G = 5  # 天气好，施工收益
loss_B = -1  # 天气不好，罚款
loss_no_work = -0.2  # 不施工窝工费
cost_of_report = 0.08  # 资料费

# 不进行预报的期望收益
expected_no_report = np.dot(P, [gain_G, loss_B])

# 进行预报后的期望收益
expected_with_report_good = P_G_given_A * gain_G + P_B_given_A * loss_B - cost_of_report
expected_with_report_bad = P_G_given_N * gain_G + P_B_given_N * loss_B - cost_of_report

expected_with_report = P_A * max(expected_with_report_good, loss_no_work - cost_of_report) + \
                       P_N * max(expected_with_report_bad, loss_no_work - cost_of_report)

# 输出概率矩阵和期望收益
print("初始概率矩阵:")
print("P(G): {:.2f}, P(B): {:.2f}".format(P[0], P[1]))
print("\n条件概率矩阵:")
print("P(A|G): {:.2f}, P(A|B): {:.2f}".format(P_A_given[0], P_A_given[1]))
print("P(N|G): {:.2f}, P(N|B): {:.2f}".format(P_N_given[0], P_N_given[1]))
print("\n后验概率矩阵:")
print("P(G|A): {:.2f}, P(B|A): {:.2f}".format(P_G_given_A, P_B_given_A))
print("P(G|N): {:.2f}, P(B|N): {:.2f}".format(P_G_given_N, P_B_given_N))
print("\n全概率:")
print("P(A): {:.2f}".format(P_A))
print("P(N): {:.2f}".format(P_N))
print("\n期望收益:")
print("不进行预报的期望收益: {:.2f}".format(expected_no_report))
print("进行预报后的期望收益: {:.2f}".format(expected_with_report))

# 比较期望收益，选择最优决策
if expected_with_report > expected_no_report:
    print("最佳决策是: 进行预报")
else:
    print("最佳决策是: 不进行预报")

初始概率矩阵:
P(G): 0.30, P(B): 0.70

条件概率矩阵:
P(A|G): 0.80, P(A|B): 0.10
P(N|G): 0.20, P(N|B): 0.90

后验概率矩阵:
P(G|A): 0.77, P(B|A): 0.23
P(G|N): 0.09, P(B|N): 0.91

全概率:
P(A): 0.31
P(N): 0.69

期望收益:
不进行预报的期望收益: 0.80
进行预报后的期望收益: 0.91
最佳决策是: 进行预报

练习3

某地区有甲、乙、丙三家食品厂生产同一种食品，有一千个用户（或购货点），假定在研究期间无新用户加入也无老用户退出，只有用户的转移，已知 2006 年 5 月份有 500 户是甲厂的顾客；400 户是乙厂的顾客；100 户是丙厂的顾客。6 月份，甲厂有 400 户原来的顾客，上月的顾客有 50 户转乙厂，50 户转丙厂；乙厂有 300 户原来的顾客，上月的顾客有 20 户转甲厂，80 户转丙厂；丙厂有 80 户原来的顾客，上月的顾客有 10 户转甲厂，10 户转乙厂。计算其状态转移概率与平稳分布。
由题意得 6 月份顾客转移表，如下：

厂名	甲	乙	丙	合计
甲	400	50	50	500
乙	20	300	80	400
丙	10	10	80	100
合计	430	360	210	1000

import numpy as np

# 原始数据矩阵
data = np.array([
    [400, 50, 50],
    [20, 300, 80],
    [10, 10, 80]
])

# 每行的合计数
row_totals = np.array([500, 400, 100]).reshape(-1, 1)

# 将数据矩阵每行除以合计数，归一化处理为转移矩阵
P = data / row_totals

print("转移矩阵 P:")
print(np.round(P, 2))

# 计算平稳分布
eigvals, eigvecs = np.linalg.eig(P.T)
eigvec1 = eigvecs[:, np.isclose(eigvals, 1)]

# 将特征向量归一化以使其总和为1
stationary = eigvec1 / eigvec1.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)

print("\n平稳分布:")
print(stationary)
#[0.29 0.29 0.43]

练习4

计算机运行状态的预测。设一随机系统状态空间, 总共有4种状态, 记录观测系统所处状态如下:

\begin{array}{llllllllll} 4 & 3 & 2 & 1 & 4 & 3 & 1 & 1 & 2 & 3 \\ 2 & 1 & 2 & 3 & 4 & 4 & 3 & 3 & 1 & 1 \\ 1 & 3 & 3 & 2 & 1 & 2 & 2 & 2 & 4 & 4 \\ 2 & 3 & 2 & 3 & 1 & 1 & 2 & 4 & 3 & 1 \end{array}

$\begin{array}{llllllllll} 4 & 3 & 2 & 1 & 4 & 3 & 1 & 1 & 2 & 3 \\ 2 & 1 & 2 & 3 & 4 & 4 & 3 & 3 & 1 & 1 \\ 1 & 3 & 3 & 2 & 1 & 2 & 2 & 2 & 4 & 4 \\ 2 & 3 & 2 & 3 & 1 & 1 & 2 & 4 & 3 & 1 \end{array}$

若该系统可用马尔科夫模型描述, 且当前系统处于状态4 , 则下一步最可能是状态为何？

求解思路：

求解一步转移概率: 从状态 $i$ 到状态 $j$ 的概率的估计值为

{\hat{p}}_{i j} = \frac{n_{i j}}{\sum_{j = 1}^{4} n_{i j}}

$\hat{p}_{i j}=\frac{n_{i j}}{\sum_{j=1}^4 n_{i j}}$

例如, 从状态1可以转移到状态 $1 、 2 、 3 、 4$ , 那么状态1转移到其他状态的总次数, 就是观测到的数据中, 出现 “ 11 ”" 12 ”"13" "14" 的总次数, 也就是公式中的分母。
从状态 $i$ 到状态 $j$ 的次数统计（不需要手算, 代码部分会算）

	1	2	3	4	$\sum_{j=1}^4 n_{ij}$
1	4	4	1	1	10
2	3	2	4	2	11
3	4	4	2	1	11
4	0	1	4	2	7

求出一步转移矩阵:

P = [\begin{array}{cccc} 2 / 5 & 2 / 5 & 1 / 10 & 1 / 10 \\ 3 / 11 & 2 / 11 & 4 / 11 & 2 / 11 \\ 4 / 11 & 4 / 11 & 2 / 11 & 1 / 11 \\ 0 & 1 / 7 & 4 / 7 & 2 / 7 \end{array}]

$P=\left[\begin{array}{cccc} 2 / 5 & 2 / 5 & 1 / 10 & 1 / 10 \\ 3 / 11 & 2 / 11 & 4 / 11 & 2 / 11 \\ 4 / 11 & 4 / 11 & 2 / 11 & 1 / 11 \\ 0 & 1 / 7 & 4 / 7 & 2 / 7 \end{array}\right]$

其中第 $i$ 行第 $j$ 列的数, 代表着从状态 $i$ 一步转移到状态 $j$ 的概率的估计值。题目说当前状态为 4 , 则下一步的状态最有可能是3。

则第 $\mathrm{n}$ 步的状态的概率分布为:
初始概率分布行向量（即初始时每种状态出现的概率）记做 $P^{(0)}$ , 一步转移概率矩阵记做 $P$ 。

P^{(n)} = P^{(0)} P^{n}

$P^{(n)}=P^{(0)} P^n$

注意是把矩阵 $P$ 乘 $n$ 次。
系统初始化后运行到第2步的概率分布: $P^{(2)}=P^{(0)} P^2=[0.291,0.289,0.271,0.149]$
系统初始化后运行到第5步的概率分布: $P^{(5)}=P^{(0)} P^5=[0.294,0.289,0.268,0.149]$

代码求解
求一步转移矩阵:一般需要根据已知数据，求出一步转移矩阵的估计值；
需要遍历已知数据，统计每一组相邻两个状态的数量（找出有多少个 "1 1"，“12",..... )；
如果题目所定的系统状态是文字（例如张三的四种状态），那么在写代码时，可以自行定义不同状态对应的数字，给出一步转移概率和转移矩阵:

{\hat{p}}_{i j} = \frac{n_{i j}}{\sum_{j = 1}^{4} n_{i j}}

$\hat{p}_{i j}=\frac{n_{i j}}{\sum_{j=1}^4 n_{i j}}$

import numpy as np

# 原始数据矩阵
a0 = np.array([
    [4, 3, 2, 1, 4, 3, 1, 1, 2, 3],
    [2, 1, 2, 3, 4, 4, 3, 3, 1, 1],
    [1, 3, 3, 2, 1, 2, 2, 2, 4, 4],
    [2, 3, 2, 3, 1, 1, 2, 4, 3, 1]
])

# 将矩阵展平为一维数组
a = a0.flatten()

# 初始化频数矩阵
freq_matrix = np.zeros((4, 4), dtype=int)

# 统计频数
for k in range(len(a) - 1):
    i, j = a[k], a[k + 1]
    freq_matrix[i - 1, j - 1] += 1

# 打印频数矩阵
print("频数矩阵:")
print(freq_matrix)

# 计算状态转移矩阵 P
sums = freq_matrix.sum(axis=1)
P = freq_matrix / sums[:, None]

# 打印状态转移矩阵 P
print("\n状态转移矩阵P:")
print(np.round(P, 2))

# 计算 P 的 2 次方和 5 次方
P2 = np.linalg.matrix_power(P, 2)
P5 = np.linalg.matrix_power(P, 5)

# 打印 P 的 2 次方和 5 次方
print("\nP 的 2 次方:")
print(np.round(P2, 2))
print("\nP 的 5 次方:")
print(np.round(P5, 2))

# 计算平稳分布
eigvals, eigvecs = np.linalg.eig(P.T)
stationary = eigvecs[:, np.isclose(eigvals, 1)]
stationary = stationary / stationary.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)

# 打印平稳分布
print("\n平稳分布:")
print(stationary)

P 的 5 次方:
[[0.29 0.29 0.27 0.15]
 [0.29 0.29 0.27 0.15]
 [0.29 0.29 0.27 0.15]
 [0.29 0.29 0.27 0.15]]

平稳分布:
[0.29 0.29 0.27 0.15]

练习5

考虑中国移动，中国联通和中国电信三家运营商，他们的用户可以互相携号转网，已经当前3家运营商的市场份额，而且也能测试出用户转网的可能性，那么将来3家运营商的市场份额情况如何，即利用当前已知的两项数据，分别是当前的市场份额、用户接下来使用运营商的可能性（即转移概率矩阵），则可预测将来3家运营商的市场份额情况。当前3家运营商分别的市场占比为0.55,0.25和0.2，但是有携号转网政策后，用户可自由的切换运营商，当前从调查数据中可得到，移动用户预期下年还会使用移动的比例是80%，使用电信的比例是15%，联通的比例是5%。电信用户接下来会使用移动的比例是20%，继续使用电信的比例是78%，使用联通的比例是2%。联通用户接下来使用移动的比例是5%，电信为20%，继续使用联通的比例是75%。现希望预期后续比如10年后，3家运营商的市场份额情况为何。


import numpy as np

# 初始市场份额
initial_distribution = np.array([0.55, 0.25, 0.2])

# 状态转移矩阵
P = np.array([
    [0.80, 0.15, 0.05],
    [0.20, 0.78, 0.02],
    [0.05, 0.20, 0.75]
])

# 计算10年后的市场份额分布
distribution_10_years = initial_distribution.dot(np.linalg.matrix_power(P, 10))

# 状态名
states = ["移动", "电信", "联通"]

print("10年后的市场份额分布:")
for state, share in zip(states, np.round(distribution_10_years, 2)):
    print(f"{state}: {share}")

# 计算平稳分布
eigvals, eigvecs = np.linalg.eig(P.T)
eigvec1 = eigvecs[:, np.isclose(eigvals, 1)]

# 将特征向量归一化以使其总和为1
stationary = eigvec1 / eigvec1.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)

print("\n平稳分布:")
for state, share in zip(states, stationary):
    print(f"{state}: {share}")

10年后的市场份额分布:
移动: 0.45
电信: 0.42
联通: 0.13

平稳分布:
移动: 0.45
电信: 0.42
联通: 0.12

posted @ 2024-06-16 12:58 郝hai 阅读(25) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 运筹学练习Python精解——排队论

· 决策论——贝叶斯决策的R和Python实现（二）

· 精通金融模式识别-全-

· 贝叶斯优化实战（全）

· UCB Data100：数据科学的原理和技巧：第二十一章到第二十六章

阅读排行：
· 全程不用写代码，我用AI程序员写了一个飞机大战
· DeepSeek 开源周回顾「GitHub 热点速览」
· 记一次.NET内存居高不下排查解决与启示
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布：重大改进与新特性概览！

公告

昵称：郝hai
园龄： 2年10个月
粉丝： 96
关注： 0

+加关注

2025年3月

日

一

二

三

四

五

六

haohai9309

格物致知，知行合一！

运筹学练习Python精解——决策论

练习1

练习2

练习3

练习4

练习5

公告

搜索

常用链接

随笔分类

随笔档案

阅读排行榜

评论排行榜

推荐排行榜

最新评论