python - 如何对方程进行蒙特卡罗分析？

Question

给定一个依赖于多个变量的函数，每个变量都有一定的概率分布，我如何进行蒙特卡罗分析以获得函数的概率分布。理想情况下，随着参数数量或迭代次数的增加，我希望解决方案具有高性能。

例如，我为 total_time 提供了一个方程式，它取决于许多其他参数。

import numpy as np
import matplotlib.pyplot as plt

size = 1000

gym = [30, 30, 35, 35, 35, 35, 35, 35, 40, 40, 40, 45, 45]

left = 5
right = 10
mode = 9
shower = np.random.triangular(left, mode, right, size)

argument = np.random.choice([0, 45], size, p=[0.9, 0.1])

mu = 15
sigma = 5 / 3
dinner = np.random.normal(mu, sigma, size)

mu = 45
sigma = 15/3
work = np.random.normal(mu, sigma, size)

brush_my_teeth = 2

variables = gym, shower, dinner, argument, work, brush_my_teeth
for variable in variables:
    plt.figure()
    plt.hist(variable)
plt.show()


def total_time(variables):
    return np.sum(variables)

健身房

淋浴

晚餐

争论

工作

刷牙

Answer 1

现有答案的想法是正确的，但我怀疑您是否想像 nicogen 那样对 size 中的所有值求和。

我假设您选择了一个相对较大的 size 来展示直方图中的形状，而您想要从每个类别中求和一个值。例如，我们要计算每个事件的一个实例的总和，而不是 1000 个实例。

第一个代码块假定您知道您的函数是求和，因此可以使用快速 numpy 求和来计算求和。

import numpy as np
import matplotlib.pyplot as plt

mc_trials = 10000

gym = np.random.choice([30, 30, 35, 35, 35, 35, 
                    35, 35, 40, 40, 40, 45, 45], mc_trials)
brush_my_teeth = np.random.choice([2], mc_trials)
argument = np.random.choice([0, 45], size=mc_trials, p=[0.9, 0.1])
dinner = np.random.normal(15, 5/3, size=mc_trials)
work = np.random.normal(45, 15/3, size=mc_trials)
shower = np.random.triangular(left=5, mode=9, right=10, size=mc_trials)

col_per_trial = np.vstack([gym, brush_my_teeth, argument,
           dinner, work, shower])

mc_function_trials = np.sum(col_per_trial,axis=0)

plt.figure()
plt.hist(mc_function_trials,30)
plt.xlim([0,200])
plt.show()

如果您不了解您的函数，或者无法轻松地将其重铸为 numpy 逐元素矩阵运算，您仍然可以像这样循环:

def total_time(variables):
        return np.sum(variables)

mc_function_trials = [total_time(col) for col in col_per_trial.T]

您询问有关获取“概率分布”的问题。像我们上面所做的那样获取直方图并不完全适合你。它为您提供了视觉表示，但不是分布函数。为了得到这个函数，我们需要使用核密度估计。 scikit-learn 有一个 jar 头 function and example就是这样做的。

from sklearn.neighbors import KernelDensity
mc_function_trials = np.array(mc_function_trials)
kde = (KernelDensity(kernel='gaussian', bandwidth=2)
       .fit(mc_function_trials[:, np.newaxis]))

density_function = lambda x: np.exp(kde.score_samples(x))

time_values = np.arange(200)[:, np.newaxis]
plt.plot(time_values, density_function(time_values))

现在您可以计算总和小于 100 的概率，例如:

import scipy.integrate as integrate
probability, accuracy = integrate.quad(density_function, 0, 100)
print(probability)
# prints 0.15809