侧边栏壁纸
博主头像
LittleAO的学习小站 博主等级

在知识的沙漠寻找绿洲

  • 累计撰写 125 篇文章
  • 累计创建 27 个标签
  • 累计收到 0 条评论

目 录CONTENT

文章目录
技术分享

数学建模系列(6)——概率模型

LittleAO
LittleAO
2023-08-22 / 0 评论 / 0 点赞 / 10 阅读 / 0 字
温馨提示：
本文最后更新于2023-11-10，若内容或图片失效，请留言反馈。 部分素材来自网络，若不小心影响到您的利益，请联系我们删除。

第六章 概率模型

这方面需要大量的概率论的知识，详情请见：

概率论与数理统计 | 标签 | LittleAO学习小站

决策模型

什么是决策？

决策者为了实现预定的目标，根据一定的条件，提出实现目标的各种行动方案，并针对每一方案在实施过程中可能面临的客观状态，运用适当的决策准则与方法，比较各方案的优劣，从中选出最优或较满意的方案加以实施的完整过程。

决策问题主要分三类：

  • 确定性决策：决策环境和选择结果都是完全确定的。
  • 不确定性决策：决策之将对发生结果的概率一无所知，只能凭决策者的主观倾向进行决策。
  • 风险性决策：决策的环境不是完全确定的，发生的概率是已知的。

在平常的建模过程，确定性建模不需要过多赘述。主要放在不确定性决策和风险性决策上面。

例题 工厂销售策略（不确定性决策）

设某工厂是按批生产并按批销售某种产品，每件产品的成本为30元，批发价格为35元。若每月生产的产品当月销售不完，则每件损失1元。工厂每批生产10件，最大月生产能力是40件，决策者可选择的生产方案为0,10,20,30,40五种。假设决策者不清楚其产品的需求情况，试问决策者应如何决策？

在这道例题中，决策者可供选择的方案有五种，即生产的产品0,10,20,30,40件。将这些决策记为S_i，为策略集合。

对应的也只有5种销售情况，为0,10,20,30,40，但不知道他们发生的概率，所以这个决策是不确定性决策。将这些时间记为记为E_j，为事件集合。

每一个策略-事件对都可以计算出相应的收益值或损失值。将每个决策和事件进行一一匹配，那么我们可以很清楚的得出一个表格：

策略\事件 0 10 20 30 40
0 0 0 0 0 0
10 -10 50 50 50 50
20 -20 40 100 100 100
30 -30 30 90 150 150
40 -40 20 80 140 200

称该表为收益矩阵。

对于做出的决策有以下五个分类：

  • 悲观决策准则：亦称为保守决策准则，考虑在最坏的情况下来争取最好的可能。

  • 乐观决策准则：选择好中最好的决策。

  • 等可能性准则：因为没有概率作为依据，那么决策者只能认为每一种情况发生的机会是均等的。决策者会计算收益的期望值，然后选取已它对应的策略作为决策策略。

  • 最小机会损失决策准则：亦称为最小遗憾值决策准则。具体做法就是将收益矩阵中的各元素变换为机会损失值。其含义是当某一事件发生后，由于决策者没有选用收益最大的策略而形成的损失值。

  • 现实主义决策标准：将悲观决策准则和乐观决策准则进行总和，取一个乐观系数\alpha，用以下关系是来表示：

    H_i=\alpha a_i\max+(1-\alpha)a_i\min

    其中a_i\max和a_i\min分别表示第i个策略可能得到的最大收益值和最小收益值。然后选择最大的H_i值为最终决策。

用悲观决策准则来进行决策，也就是选取：

\max\{0,-10,-20,-30,-40\}=0

对应的决策，也就是一件也不生产。解释起来就是只要不生产，就不会赔钱。

用乐观决策准则来进行决策，也就是选取：

\max\{0,50,100,150,200\}=200

对应的决策也就是生产40件，解释起来是生产的越多，我们获利的利润就越大。

用等可能性准则来看，我们需要计算每个决策的期望，然后选取期望最大的那一个：

\max\{E(S_i)\} = S_5=80

对应的决策就是生产50件。

使用最小机会损失准则，算出的损失矩阵如图：

策略\事件 0 10 20 30 40
0 0 50 100 150 200
10 10 0 50 100 150
20 20 10 0 50 100
30 30 20 10 0 50
40 40 30 20 10 0 
import numpy as np
# 根据收益矩阵计算损失矩阵

def loss_matrix(profit_matrix):
    """根据收益矩阵计算损失矩阵

    Parameters
    ----------
    profit_matrix : np.array
        收益矩阵

    Returns
    -------
    np.array
        损失矩阵
    """
    loss_matrix = np.zeros(profit_matrix.shape)
    # 找出每一行的最大值
    max_value = np.max(profit_matrix, axis=1)
    # 更新损失矩阵，每一列为每一行的最大值
    for i in range(profit_matrix.shape[0]):
        loss_matrix[:, i] = max_value[i]
    # 损失矩阵的每个元素减去收益矩阵的每个元素
    loss_matrix = loss_matrix - profit_matrix
    return loss_matrix

def confirm(loss_matrix):
    """使用最小机会损失准则确定最优策略

    Parameters
    ----------
    loss_matrix : np.array
        损失矩阵

    return
    ------
    int 
        最优策略编号
    """
    # 找出每一行的最大值
    min_value = np.max(loss_matrix, axis=1)
    # 找出最大值中的最小值
    min_value = np.min(min_value)
    # 找出最小值所在的行
    min_index = np.where(loss_matrix == min_value)[0][0]
    return min_index + 1

if __name__ == '__main__':
    profit_matrix = np.array(
        [[0, 0, 0, 0, 0],
         [-10, 50, 50, 50, 50],
         [-20, 40, 100, 100, 100],
         [-30, 30, 90, 150, 150],
         [-40, 20, 80, 140, 200]])
    loss_matrix = loss_matrix(profit_matrix)
    # print(loss_matrix)
    print("选择编号为{}的策略".format(confirm(loss_matrix)))

选取策略的方法：

\min\{200, 150,100,50,40\} = 40

意为选择策略5的损失最小。

使用现实主义决策标准，我们首先假设\alpha = 1/3，然后计算H_i值：

import numpy as np
# 现代主义决策准则

def realistic(profit_matrix, alpha=0.5):
    """现代主义决策准则

    Parameters
    ----------
    profit_matrix : np.array
        收益矩阵
    alpha : float, optional
        期望收益权重, 默认0.5

    Returns
    -------
    int
        最优策略编号
    """
    # 计算各决策的H值
    max_profit = np.max(profit_matrix, axis=1)
    min_profit = np.min(profit_matrix, axis=1)
    H = alpha * max_profit + (1 - alpha) * min_profit
    # print(H)
    # 返回最优策略编号
    return np.argmax(H) + 1

if __name__ == '__main__':
    profit_matrix = np.array(
        [[0, 0, 0, 0, 0],
         [-10, 50, 50, 50, 50],
         [-20, 40, 100, 100, 100],
         [-30, 30, 90, 150, 150],
         [-40, 20, 80, 140, 200]])
    print(realistic(profit_matrix, 1/3))

经过计算：

H=\{0,10,20,30,40\}

故选择决策5作为最终决策。

在进行不确定性决策时，要因人、因地、因时的选择决策模型。

例题 化工厂创新（风险决策）

风险决策模型相比于不确定性决策模型多了不同事件发生的概率。风险决策准则分为三大类：

  • 最大期望收益决策准则：将计算每个策略——事件对的期望收益值，然后选取期望收益最大的决策。
  • 最小机会损失决策模型：计算每个策略——事件对的机会损失值，选取损失值最小的决策。
  • 决策树：将每个可能的选项和收益化成树状图进行分析。决策树分为三个要素：决策点（用方框表显示）、状态点（用圆圈表示）、树梢（用三角表示，代表损益值）。

下面的例题来详细解释决策树的用法：

某一化工原料厂，由于某项工艺不甚好，产品成本高。在价格中等水平时无利可图，在价格低落时要亏本，只有在价格高时才赢利，且赢利也不多。

现企业考虑进行技术革新，取得新工艺的途径有两种，一是自行研究，研究成功的可能是0.6，二是购买专利，估计购买谈判成功的可能性是0.8。

不论是研究成功还是谈判成功，生产规模有两种考虑方案，一是产量不变，二是产量增加。若研究失败或者谈判失败，则仍然采用原工艺进行生产，生产保持不变。

根据市场预测，今后五年内产品跌价的可能性是0.1，保持中等水平的可能性是0.5，涨价的可能性是0.4。决策收益值见下表。现在企业需要决定：购买专利还是自行研究。

方案 按原工艺生产 购买专利成功(0.8) 购买专利成功(0.8) 自行研究成功(0.6) 自行研究成功(0.6)
概率 产量不变 增加产量 产量不变 增加产 量
价格低落 (0.1) -1000 -2000 -3000 -2000 -3000
中等(0.5) 0 500 500 0 -2500
高涨(0.4) 1000 1500 2500 2000 6000

我们首先画出决策树：

graph LR
  1 -->|购买专利| 2((2))
  1 -->|自行研究| 3((3))
	2 --->|失败 0.2| 4((4))
	2 -->|成功 0.8| 5[5]
	3 -->|成功 0.6| 6[6]
	3 --->|失败 0.4| 7((7))
	5 --> |产量不变|8((8))
	5 --> |增加产量|9((9))
	6 --> |产量不变|10((10))
	6 --> |增加产量|11((11))
	4 --> |0.1|12[/-1000\]
	4 --> |0.5|13[/0\]
	4 --> |0.4|14[/1000\]
	8 --> |0.1|15[/-2000\]
	8 --> |0.5|16[/500\]
	8 --> |0.4|17[/1500\]
	9 --> |0.1|18[/-3000\]
	9 --> |0.5|19[/500\]
	9 --> |0.4|20[/2500\]
	10 --> |0.1|21[/-2000\]
	10 --> |0.5|22[/0\]
	10 --> |0.4|23[/2000\]
	11 --> |0.1|24[/-3000\]
	11 --> |0.5|25[/-2500\]
	11 --> |0.4|26[/6000\]
	7 --> |0.1|27[/-1000\]
	7 --> |0.5|28[/0\]
	7 --> |0.4|29[/1000\]

紧接着，我们通过从左往右算，算出每个状态点的数学期望，当算到决策点时，我们选择数学期望最优的决策，将其他决策删掉：

graph LR
  1 -->|购买专利| 2((2))
  1 -->|自行研究| 3((3))
	2 --->|失败 0.2| 4((4))
	2 -->|成功 0.8| 5[5]
	3 -->|成功 0.6| 6[6]
	3 --->|失败 0.4| 7((7))
	5 --> |产量不变|8((8))
	5 --> |增加产量|9((9))
	6 --> |产量不变|10((10))
	6 --> |增加产量|11((11))
	subgraph a1[300]
	4 --> |0.1|12[/-1000\]
	4 --> |0.5|13[/0\]
	4 --> |0.4|14[/1000\]
	end
	subgraph a2[650]
	8 --> |0.1|15[/-2000\]
	8 --> |0.5|16[/500\]
	8 --> |0.4|17[/1500\]
	end
	subgraph a3[900]
	9 --> |0.1|18[/-3000\]
	9 --> |0.5|19[/500\]
	9 --> |0.4|20[/2500\]
	end
	subgraph a4[600]
	10 --> |0.1|21[/-2000\]
	10 --> |0.5|22[/0\]
	10 --> |0.4|23[/2000\]
	end
	subgraph a5[850]
	11 --> |0.1|24[/-3000\]
	11 --> |0.5|25[/-2500\]
	11 --> |0.4|26[/6000\]
	end
	subgraph a6[300]
	7 --> |0.1|27[/-1000\]
	7 --> |0.5|28[/0\]
	7 --> |0.4|29[/1000\]
	end

删掉不必要的状态：

graph LR
  1 -->|购买专利| 2((2))
  1 -->|自行研究| 3((3))
	2 --->|失败 0.2| 4((4))
	2 -->|成功 0.8| 5[5]
	3 -->|成功 0.6| 6[6]
	3 --->|失败 0.4| 7((7))
	subgraph a1[300]
	4 --> |0.1|12[/-1000\]
	4 --> |0.5|13[/0\]
	4 --> |0.4|14[/1000\]
	end
	subgraph a3[900]
	5 --> |增加产量|9((9))
	9 --> |0.1|18[/-3000\]
	9 --> |0.5|19[/500\]
	9 --> |0.4|20[/2500\]
	end
	subgraph a5[850]
	6 --> |增加产量|11((11))
	11 --> |0.1|24[/-3000\]
	11 --> |0.5|25[/-2500\]
	11 --> |0.4|26[/6000\]
	end
	subgraph a6[300]
	7 --> |0.1|27[/-1000\]
	7 --> |0.5|28[/0\]
	7 --> |0.4|29[/1000\]
	end

然后接着计算数学期望：

graph LR
  1 -->|购买专利| 2((2))
  1 -->|自行研究| 3((3))
	
	subgraph b1[780]
	2 --->|失败 0.2| 4((4))
	2 -->|成功 0.8| 5[5]
	subgraph a1[300]
	4 --> |0.1|12[/-1000\]
	4 --> |0.5|13[/0\]
	4 --> |0.4|14[/1000\]
	end
	subgraph a3[900]
	5 --> |增加产量|9((9))
	9 --> |0.1|18[/-3000\]
	9 --> |0.5|19[/500\]
	9 --> |0.4|20[/2500\]
	end
	end
	subgraph b2[630]
	3 -->|成功 0.6| 6[6]
	3 --->|失败 0.4| 7((7))
	subgraph a5[850]
	6 --> |增加产量|11((11))
	11 --> |0.1|24[/-3000\]
	11 --> |0.5|25[/-2500\]
	11 --> |0.4|26[/6000\]
	end
	subgraph a6[300]
	7 --> |0.1|27[/-1000\]
	7 --> |0.5|28[/0\]
	7 --> |0.4|29[/1000\]
	end
	end

删掉不必要的状态：

graph LR
  1 -->|购买专利| 2((2))
	
	subgraph b1[780]
	2 --->|失败 0.2| 4((4))
	2 -->|成功 0.8| 5[5]
	subgraph a1[300]
	4 --> |0.1|12[/-1000\]
	4 --> |0.5|13[/0\]
	4 --> |0.4|14[/1000\]
	end
	subgraph a3[900]
	5 --> |增加产量|9((9))
	9 --> |0.1|18[/-3000\]
	9 --> |0.5|19[/500\]
	9 --> |0.4|20[/2500\]
	end
	end

求解得到最终的决策应该为购买专利，且增加产量。

例题 购买飞机（多属性决策）

某航空公司在国际市场上购买飞机，按6个决策指标对不同型号的飞机进行综合评价，这6个指标是：最大速度、最大范围、最大负载、价格、可靠性、灵敏度。现在4种型号的飞机可供选择，具体指标值见表。

最大速度 (马赫) 最大范围 (千米) 最大负载 (千克) 费用 (百万元) 可靠性 灵敏度
1 2.0 1500 20000 5.5 一般 很高
2 2.5 2700 18000 6.5 低 一般
3 1.8 2000 21000 4.5 高 高
4 2.2 1800 20000 5.0 一般 一般

多属性决策构建决策矩阵，我们将可靠性和灵敏度量化。得到的决策矩阵为：

X=(x_{y})_{4\times6}=\begin{bmatrix}2.0&1500&20000&5.5&5&9\\2.5&2700&18000&6.5&3&5\\1.8&2000&21000&4.5&7&7\\2.2&1800&20000&5.0&5&5\end{bmatrix}

接着我们需要将不同的指标进行标准化，在这里我们用到的方法是向量归一化法：

在决策矩阵X=(x_{ij})_{m\times n}中，令：

y_{ij}=\frac{{x_{ij}}}{\sqrt{\sum\limits_{i=1}^{m}x^2{}_{ij}}}\left(1\leq i\leq m,1\leq j\leq n\right)

所得到的矩阵Y，称为向量归一化矩阵，每一个指标值均满足0\le y_{ij}\le1，并且正、逆向指标的方向没有发生变化。

Y= \begin{bmatrix} 0.4671& 0.3662& 0.5056& 0.5069& 0.4811& 0.6708\\ 0.5839& 0.6591& 0.4550 & 0.5990 & 0.2887& 0.3727\\ 0.4204& 0.4882& 0.5308& 0.4147& 0.6736& 0.5217\\ 0.5139& 0.4394& 0.5056& 0.4608& 0.4811& 0.3727 \end{bmatrix} 

# 向量归一化
import numpy as np

def normalize(v):
    """归一化向量，将每一列的值的开方除以该列的每个元素的平方和的开方

    Parameters
    ----------
    v : numpy.array
        需要归一化的向量

    Returns
    -------
    numpy.array
        归一化后的向量
    """
    return v / np.sqrt(np.sum(v ** 2, axis=0))

if __name__ == '__main__':
    a = np.array([[2, 1500, 20000, 5.5, 5, 9],
                  [2.5, 2700, 18000, 6.5, 3, 5],
                  [1.8, 2000, 21000, 4.5, 7, 7],
                  [2.2, 1800, 20000, 5.0, 5, 5]])

或者使用极差变换法：

Y=\begin{bmatrix}0.28&0&0.67&0.50&0.51&1.00\\1.00&1.00&0&0&0&0\\0&0.42&1.00&1.00&1.00&0.50\\0.57&0.52&0.67&0.25&0.50&0\end{bmatrix} 

def normalize(v):
    """极差变换法进行归一化

    Parameters
    ----------
    v : np.array
        输入的向量

    return
    ------
    np.array
        归一化后的向量
    """
    # 求极差
    r = np.max(v, axis=0) - np.min(v, axis=0)
    # 极差变换
    v = (v - np.min(v, axis=0)) / r
    return v

或者使用线性比例法：

def normalize(v):
    """使用线性比例法进行向量归一化

    Parameters
    ----------
    v : np.array
        输入的向量

    return
    ------
    np.array
        归一化后的向量
    """
    max_v = np.max(v, axis=0)
    return v / max_v

Y=\begin{bmatrix}0.80&0.56&0.95&0.82&0.71&1.00\\1.00&1.00&0.86&0.69&0.43&0.56\\0.72&0.74&1.00&1.00&1.00&0.78\\0.88&0.67&0.95&0.90&0.71&0.56\end{bmatrix}

得到了决策矩阵，做出决策我们可以采用下面两种方法：

  • 线性加权：对不同决策的指标加权，然后测出各方案的综合权值，来选取可行方案的排序。我们可以使用层次分析法来确定不同的权值。

  • TOPSIS法：以靠近理想解和远离负理想解两个基准作为评价各可行方案的依据。每列最大为理想解，每列最小为负理想解。测度方法如下：

    S_i^*=\sqrt{\sum_{j=1}^n\left(x_{ij}-x_j^*\right)^2}\quad S_i^-=\sqrt{\sum_{j=1}^n\left(x_{ij}-x_j^-\right)^2}

    相对贴近度：C_i^*=\frac{S_i^-}{S_i^-+S_i^*}

我们使用理想解法和向量归一化法来进行决策，假设权值W=(0.2,0.1,0.1,0.1,0.2,0.3)^T，首先我们算出加权标准化矩阵，将权值和对应的列相乘，然后用理想解法得出相对贴近度为：

C^*=(0.6433,0.2871,0.5972,0.3016)

由此可见，决策1最优。

def ideal_solution(v, w):
    """使用理想解法做出决策

    Parameters
    ----------
    v : np.array
        输入的向量
    w : np.array
        权重
    return
    ------
    int 
        最优解的索引
    """
    # 求理想解，为每一列的最大值
    v = v * w.reshape(1, -1)
    print(v)
    ideal = np.max(v, axis=0)
    print(ideal)
    # 求负理想解，为每一列的最小值
    negative_ideal = np.min(v, axis=0)
    print(negative_ideal)
    # 求每个向量到理想解的距离
    distance_ideal = np.sqrt(np.sum((v - ideal) ** 2, axis=1))
    print(distance_ideal)
    # 求每个向量到负理想解的距离
    distance_negative_ideal = np.sqrt(
        np.sum((v - negative_ideal) ** 2, axis=1))
    # 求每个向量到理想解的距离与到负理想解的距离的比值
    ratio = distance_negative_ideal / \
        (distance_ideal + distance_negative_ideal)
    print(ratio)
    # 返回最大值的索引
    return np.argmax(ratio) + 1

随机存储模型

在之前的学习中，我们已经接触过了确定性存储模型：

数学建模系列(2)——初等模型 | LittleAO学习小站

下面我们来探究随机性存储模型

例题 报童问题

上面的例子是一个连续且确定的模型，而下面要讨论的模型是离散且随机的模型，下面看这样一道例题：

报童每天清晨用每份2元的价格从报社买进一批报纸后，在报亭以每份4元零售，到晚上将没有卖掉的报纸退回，得到每份1元的补偿。这样，他每售出一份报纸可赚2元，每退回一份要赔1元。报童每天应该买进多少报纸，当然取决于他每天能卖出多少报纸。显然，卖出的数量是随机的，因而报童每天的收入也是随机的。如果已经由统计资料或经验判断得到了报纸需求量的概率分布。那么，所谓报童的诀窍是指，他每天应该买进多少份报纸，从长期来看可以获得最高的日均利润。

由统计资料得出的需求量的概率分布如图：
需求量/单位 0 1 2 3 4 5
概率/分布 0 0.1 0.3 0.4 0.1 0.1

我们记需求量取值为r时的概率为f(r)(r=1,2,3,...,5)。假设每100份报纸记为一单位，设售出一单位报纸的获利为s_1，因剩余而退回1单位的报纸损失为s_2。设报童每天购进q单位报纸，每天的需求量r，当供不应求为r>q，其售出的获利为s_1q，没有损失。而供过于求的获利为s_1r，因剩余退回的损失为s_2(q-r)，记一天的利润为s(r,q)。可得：

s\left(r,q\right)=\begin{cases}s_1r-s_2(q-r),&r\leq q\\[2ex]s_1q,&r>q\end{cases}

记s(r,q)的期望值为E(q)，有：

E(q)=\sum_{r=0}^{n}s(r,q)f(r)=\sum_{r=0}^{q}\Big[s_{1}r-s_{2}(q-r)\Big]f(r)+\sum_{r=q+1}^{n}s_{1}qf(r)

问题归结为在f(r)的条件下，求购进报纸的最优数量q，使得日均利润E(q)最大。

使用差分的方法：

设：

\Delta E(q)=E(q+1)-E(q)

若q<r,q增加到q+1,可以多售出1单位，利润增加s1;若q≥r,q增加到q+1,因剩余而多退回1单位，利润减少s2。于是：

\Delta E(q)=s_{1}P\left(r>q\right)-s_{2}P\left(r\le q\right)

由于P(r>q)=1-P(r\le q)，代入原式可得：

E(q+1)-E(q)=s_{1}\left[1-P(r\leqslant q)\right]-s_{2}P(r\leqslant q)=s_{1}-\left(s_{1}+s_{2}\right)P(r\leqslant q)

当P(r\le q)由小于s_1/(s_1+s_2)变为大于s_1/(s_1+s_2)时E(q)达到最大。

也就是

P(r\leqslant q)\geqslant\frac{s_{1}}{s_{1}+s_{2}}

成立的最小q值。

查表得最小q值为3，即每天购进3单位报纸可使日均利润E(q)达到最大。得到E(q)=450元。

在上面我们研究的是离散型需求下的模型，下面我们来讨论一下连续型随机存储模型，还是以报童售报为例：

设由统计资料和经验判断认为连续型的报纸需求量大致服从正态分布，设其平均值是260，标准差是50，即需求量为N(260,50^2)。报童仍然每天以一份报纸2元买进、4元售出、退回得到1元补偿，即每售出1份报纸获利2元，因剩余而退回1份报纸损失1。问报童每天应该购进多少份报纸，才能获得最高的日均利润？

记报纸需求量的概率密度函数为p(r)，s_1,s_2,q的含义与之前相同，可得到日均利润为:

E(q)=\int_{-\infty}^{q}\left[s_{1}r-s_{2}(q-r)\right]p(r)\mathrm{d}r+\int_{q}^{\infty}s_{1}qp(r)\mathrm{d}r

上式对 $$q求导并等于0，设r概率密度函数为F(q)=\int^q_{-\infty}p(r)\mathrm dr，得到q的最优解满足：

F(q)=\frac{s_{1}}{s_{_1}+s_{_2}} 

from scipy import stats
# 已知N(260,50^2)的F(x)=2/3，求x
print(stats.norm.ppf(2/3,260,50))

解得：x=281.536，取x=282。

随机人口模型

随机人口模型

马尔可夫链

马尔可夫链的原理有关内容我们可以从信息论与编码中找到：

信息论和编码第二章笔记 | LittleAO学习小站

下面着重介绍马尔可夫链的例题：

例题 基因遗传

生物繁殖时，一个后代随机地继承父亲两个基因中的一个和母亲两个基因中的一个，形成它的两个基因.一般两个基因中哪一个遗传下去是等概率的，所以父母的基因类型就决定了每一后代基因类型的概率。父母基因类型有全是优种DD，全是劣种RR,一优种一混种DH(父为D,母为H或父为H,母为D)，及DR,HH,HR共6种组合，对每种组合简单的计算可以得到其后代各种基因类型的概率，如表所示。

| 后代基因类型|父母基因类型 | DD | RR | DH | DR | HH | HR |
| — | — | — | — | — | — | — |
| D | 1 | 0 | \frac 12 | 0 | \frac 14 | 0 |
| H | 0 | 0 | \frac 12 | 1 | \frac 12 | \frac 12 |
| R | 0 | 1 | 0 | 0 | \frac 14 | \frac 12 |

下面以马尔可夫链为工具讨论混种繁殖，优种繁殖和近亲繁殖3个基因遗传模型。

  • 混种繁殖

    假设在混种繁殖假设在繁殖过程中用一混种与一个个体交配，所得后代仍用混种交配，如此继续下去，称为混种繁殖。建立马尔可夫链模型描述在混种繁殖下各代具有3种基因类型的概率，并讨论稳态情况。

    我们找出所有混种的概率，简历转移概率矩阵：

    \boldsymbol{P}=\begin{bmatrix}\dfrac{1}{2}&\dfrac{1}{2}&0\\\\\dfrac{1}{4}&\dfrac{1}{2}&\dfrac{1}{4}\\\\0&\dfrac{1}{2}&\dfrac{1}{2}\end{bmatrix}

    若初始时混种和优种交配，则\boldsymbol a(0)=[1,0,0]。根据马尔可夫链的方程：

    \boldsymbol a(k+1)=\boldsymbol a(k)\boldsymbol P 

    # 马尔可夫链
import numpy as np

# 生成状态转移矩阵

def generate_transition_matrix():
    transition_matrix = np.array(
        [[0.5, 0.5, 0], [0.25, 0.5, 0.25], [0, 0.5, 0.5]])
    return transition_matrix

# 生成初始状态概率向量

def generate_initial_state_vector():
    initial_state_vector = np.array([1, 0, 0])
    return initial_state_vector

if __name__ == '__main__':
    # 生成状态转移矩阵
    transition_matrix = generate_transition_matrix()
    # 生成初始状态概率向量
    initial_state_vector = generate_initial_state_vector()
    # 生成马尔可夫链
    last_state_vector = initial_state_vector
    while True:
        current_state_vector = np.dot(last_state_vector, transition_matrix)
        if np.allclose(last_state_vector, current_state_vector, rtol=1e-10):
            break
        last_state_vector = current_state_vector
    print(current_state_vector)

    我们最终可以得到稳态转移概率：

    \boldsymbol a=[\frac 14,\frac12,\frac 14]

    也就是说，经过足够多代繁殖之后，优种、混种、劣种的比例接近于1:2:1。

    有\boldsymbol P^2 >\boldsymbol O 可得该链为正则链，最终稳态概率必定为：

    \boldsymbol a=[\frac 14,\frac12,\frac 14]

  • 优种繁殖：

    我们将混种交配改为每一代都用优种交配，得到转移概率矩阵：

    \boldsymbol{P}=\left[\begin{array}{ccc}1&0&0\\\dfrac{1}{2}&\dfrac{1}{2}&0\\0&1&0\end{array}\right]

    我们计算\boldsymbol P^n得到

    \boldsymbol{P}^n=\left[\begin{array}{ccc}1&0&0\\1&0&0\\1&0&0\end{array}\right]

    该链为吸收链，无论初始用的是是什么种，最终都将称为优种。

  • 近亲繁殖

    近亲繁殖指的是同一对父母的大量后代中，随机的选取一雄一雌进行交配，产生后代，如此继续下去。

    后代的基因型有6中可能性：DD,RR,DH,DR,HH,RR，分别计算每种基因型交配过后产生的基因型的概率，构造转移状态矩阵：

    \boldsymbol{P}=\begin{bmatrix}1&0&0&0&0&0\\0&1&0&0&0&0\\\dfrac{1}{4}&0&\dfrac{1}{2}&0&\dfrac{1}{4}&0\\0&0&0&0&1&0\\\dfrac{1}{16}&\dfrac{1}{16}&\dfrac{1}{4}&\dfrac{1}{8}&\dfrac{1}{4}&\dfrac{1}{4}\\0&\dfrac{1}{4}&0&0&\dfrac{1}{4}&\dfrac{1}{2}\end{bmatrix}

    该阵也是收敛阵，最终的后代品种都会向优种或劣种收敛。在现实生产过程中优种或劣种的性状比较差，因此在近亲繁殖4-5代后必须重新选种，以防品质的下降。

0
数学建模 Python
版权归属： LittleAO
本文链接： https://bg.littleao.top/archives/mathematical_modeling-6
许可协议： 本文使用《署名-非商业性使用-相同方式共享 4.0 国际 (CC BY-NC-SA 4.0)》协议授权

评论区