厦门大学暑期数学建模培训——随机系统作业

res = 1 - (0.3 * 0.2 * 0.1)  # 使用高中常识方法解题，利用反命题
print("被射中的概率为", res)
per1 = [0.7, 0.8, 0.9]  # 被击中概率
per2 = [0.7, 0.9, 1]  # 坠机概率
strike = [1, 0, 0, 0]  # 被击中次数概率向量，第i个位置代表了被击中i次的概率
for percent in per1: 
    temp = np.zeros_like(strike).tolist()  # 设置一个相同大小的0矩阵
    temp[0] = strike[0] * (1 - percent)  # 优先设置命中0次的概率，因为不符合公式
    for i in range(1, 4):
        temp[i] = strike[i - 1] * percent + strike[i] * (1 - percent) 
        # 命中i次的概率为命中i-1次的概率乘这一次命中概率 + 命中i次的概率乘这一次没命中的概率
    strike = temp  # 将命中概率更新
strike = sum([strike[i] * per2[i - 1] for i in range(1, 4)])  # 取出命中的概率×命中坠机的概率
print("被击毁的概率", strike)

二、正态分布（sympy.norm）

题目

一道非常基础的正态分布题

思路

关键点在于：如何获取正态分布的平均值和方差，了解一个公式

$\mu + z * \sigma = direct$

$\mu$ 为平均值， $\sigma$ 为标准差，z为z-score（ $z = \frac{p - \mu}{\sigma}$ ），direct为实际值。

通过题目中给的两个点，取得z中p值，通过norm.ppf(p)其默认的 $\mu = 0, \sigma=1$ ，获取z值

利用线性代数的方程求解，取得 $\mu$ ， $\sigma$ 的实际值

$\left\{\begin{matrix} \mu + z60 * \sigma = 60 \\ \mu + z90 * \sigma = 90 \end{matrix}\right.$ 得到 $\begin{vmatrix} 1 & z60 & 60\\ 1 & z90 & 90 \end{vmatrix}$ 求解该方程组即可

ps：np.linalg.solve无法求解无解情况和非方阵情况，使用lstsq可以求解非方阵

其结果是回归系数的值，接下来是残差平方和，矩阵的秩

代码

above_90 = 359 / 10000  # 初始化高于90分的概率
below_60 = 1151 / 10000  # 初始化低于60分概率
z60 = norm.ppf(below_60)  # ppf是逆函数，将概率转换为实际分布的位置
# 当ppf的其他两个变量未输入时默认为0和1，也就是说此时的ppf能直接转换为z值
z90 = norm.ppf(1 - above_90)  # 因为高于90分的概率p在正态分布处于右侧，其概率实际为1 - p
model = np.array([[1, z60],[1, z90]])  # 设置一个矩阵
res = np.linalg.lstsq(model, np.array([60, 90]).T, rcond=None)  # 使用lstsq最小二乘法进行求解
mean, std = res[0]  # 分别对应平均值和标准差
print("平均值：%f\n标准差：%f"%(mean, std))
model = norm(mean, std)  # 设置一个正态分布的模型，其平均值和标准差已经求解
# 与上面的ppf不同，此时调用建立好的模型ppf，实现的是概率转换为实际分布的位置
res = model.ppf(0.75) # 取得实际分布位置，也就是最低分数
print("入取最低分数线：%.2f" % res)

三、蒙特卡罗法（随机数法 + 可视化）

题目

蒙特卡罗模拟法是一种基于随机抽样来解决计算问题的方法。

思路

蒙特卡罗的一个具体例子：利用随机点在一个正方形和圆中打点，计算 $\pi$ 的近似值

积分即为计算函数与x轴形成的面积

那么我们也可以用蒙特卡罗的思想，通过随机投点得到在函数下方的点这些点占所有点的比例 * 总面积即为积分值

此时蓝色与红色共同占下了整个面积，红色点占所有点的比例，就是函数曲线下的面积占总矩形面积的比例

代码

np.random.seed(42)  # 固定随机数种子以实现结果复现
random_array = np.empty((100000, 2))  # 创建100000列2行（也就是100000个点的列表）
random_array[:, 0] = np.random.uniform(0, 1, random_array.shape[0])  # 设置随机浮点数
random_array[:, 1] = np.random.uniform(0, 1, random_array.shape[0])  # 同上
count = 0  # 计数在函数下方
for i, j in random_array:
    if 1 / (1 + i ** 2) > j:  # 符合条件的点
        count += 1
print("积分面积 × 4 = ", count / 100000 * 4)
print("pi = ", np.pi)

四、补充题目（用到销售概率矩阵）

思路

没有任何技巧，纯粹的带值计算
销售状态转移概率矩阵值得就是：畅销-畅销，畅销-滞销，滞销-滞销，滞销-畅销

这几种状态间的转移概率，并且每一行的总和为1

对于第二小题，通过概率向量获得两种概率的可能，其与基础概率统计有异曲同工之妙，但是略有不同的在于此时两个矩阵直接相乘即可

举个例子：本期上升概率 = 连续上升 × 上升概率 + 下降后上升 × 下降概率

下降概率也是同理

$\begin{vmatrix} P_{up} & P_{down} \end{vmatrix}_{1 \times 2} * \begin{vmatrix} P_{upup} & P_{updown}\\ P_{downup} & P_{downdown} \end{vmatrix} _{2 \times 2} = \begin{vmatrix} P'_{up} & P'_{down} \end{vmatrix}_{1 \times 2}$

代码

data = "112122111212112211212111"
length = len(data)  # 获取长度
upup = 0; updown = 0; downdown = 0; downup = 0; up = 0; down = 0  # 初始化
for i in range(length - 1):  # 分类4项状态转移
    this_ = data[i]  # 当期状态
    next_ = data[i + 1]  # 下一期状态
    if this_ == "1":
        up += 1  # 畅销数 + 1
        if this_ == next_:
            upup += 1  # 连续畅销 + 1
        else:
            updown += 1  # 上下 + 1
    else:
        down += 1  # 滞销数 + 1
        if this_ == next_:
            downdown += 1  # 连续滞销
        else:
            downup += 1
p = np.array([upup / up, updown / up, downup / down, downdown / down]).reshape(2, 2)# 获取状态转移概率矩阵
# 状态转移由开始状态决定，例如开始为畅销，这一行的概率分别为畅销的状态转移数 / 畅销数
print("概率矩阵为", p)  

init = np.array([1, 0])  # 因为初始状态为畅销，所以概率矩阵的开始为1畅销，0滞销
for i in range(4):
    init = np.dot(init, p)  # 将1×2的概率向量和2×2的概率矩阵相乘，可得下一个状态的概率
    print(f"第{i + 1}个季度的销售状态概率：{init}，最有可能是", "畅销" if (np.argmax(init) == 0) else "滞销", sep="")

init = [up/(up + down), down/(up + down)]  # 初始化概率向量，此时概率总和为1
up = 0; down = 0
print("初始概率向量为：", init)
for i in range(1000):
    init = np.dot(init, p)  # 同第二题思路
    if np.argmax(init) == 0:  # 获取概率更高值的索引
        up += 1 
    else:
        down += 1
print("经过1000次循环\n畅销次数为:%i\n滞销次数为:%i"%(up, down))