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

给定一个依赖于多个变量的函数，每个变量都有一个特定的概率分布，我如何做蒙特卡罗分析来获得函数的概率分布。理想情况下，我希望随着参数数量或迭代次数的增加，解决方案能够高性能地运行。
作为一个例子，我为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)

体育馆

淋浴

晚餐

论点

工作

刷牙

最佳答案

现有的答案有正确的想法，但我怀疑您是否希望像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()

如果您不知道自己的函数，或者不容易重铸，那么您仍然可以像这样循环：

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

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

你问关于获得“概率分布”。像我们上面所做的那样得到柱状图并不是很适合你。它为您提供了一个可视的表示，而不是分布函数。为了得到函数，我们需要使用核密度估计。Scikit Learn有一个罐头，可以做到这一点。

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))

function and example
现在可以计算总和小于100的概率，例如：

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

关于python - 如何对方程进行蒙特卡罗分析？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/43551666/