蒙特卡洛模型

蒙特卡罗法也称统计模拟法、统计试验法。是把概率现象作为研究对象的数值模拟方法，下面讲两个例子加深理解。

• 三门问题

  你参加一个有奖竞猜节目，现在面前有三个门，只有其中一个门后有奖，其余门后为空。你随机选中一个门后，主持人打开另外两个门中的一个门后，为空，接着，你有一个换门的机会，换还是不换？

  不考虑其他因素

  该问题已经有解，如果不换，那么主持人的操作跟你就没有关系了，所以中奖概率为\(\frac{1}{3}\)；如果换的话，你刚开始选门有三种情况，有一种是最开始选中了奖品，结果换到了其他门，而有两种可能是最开始选中了没有奖品的门，最后换到了有奖的门，那么中奖的概率就为 \(\frac{2}{3}\)。假设你现在不知道这个结论，那应该如何求解呢？代码模拟，也叫蒙特卡洛模拟，关键是量大和随机（区别于枚举）。

  模拟结果如下（代码见最后）：

  

• 蒙特卡洛求非线性规划初始值

  非线性规划详见数学规划模型

  有如下的非线性规划模型：

  \[\begin{align} &min\quad f(x) = x_1^2 + x_1x_2+x_3^2\\[0.3cm] &s.t = \begin{cases} x_1 + x_2 + x_3 \leq 100\\[0.2cm] x_1^2 + x_2^2 + x_1x_2x_3 \leq 200\\[0.2cm] x_1,x_2,x_3 \geq 0 \end{cases} \end{align} \]

  在对该非线性规划模型使用 Matlab 求解之前，需要设定 \(x_1,x_2,x_3\) 的初始值，初始值的选取对非线性规划求解非常重要，这里可以用蒙特卡洛来选取较好的初始值。

Matlab 代码

三门问题模拟代码

import random
import time

def game(result):
  random.seed(time.time())
  prize = random.randrange(1, 4) # 奖品所在门号
  first = random.randrange(1, 4) # 第一次选中的门

  # 主持人打开了一个空门

  # 不换
  if first == prize:
    result[0] += 1 # 坚持得一分
  else:
    pass # 不得奖

  # 换
  if first != prize:
    result[1] += 1 # 换得一分
  else:
    pass # 不得奖



def simulation():
  result, times = [0.0, 0.0], 1000*1000

  for i in range(times):
    game(result)

  print("共计模拟游戏 %d 次\n  换门中奖 %d 次，概率为 %lf\n  不换中奖 %d 次，概率为 %lf" 
  % (times, result[1], result[1]/times, result[0], result[0]/times))
  

if __name__ == "__main__":
    simulation()