如何对方程进行 Monte Carlo 分析？
How can I do a Monte Carlo analysis on an equation?
给定一个依赖于多个变量的函数，每个变量都有一定的概率分布，我如何进行 Monte Carlo 分析以获得函数的概率分布。理想情况下，随着参数数量或迭代次数的增加，我希望解决方案具有高性能。
例如，我为 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)
健身房
淋浴
晚餐
参数
工作
brush_my_teeth
您尝试过简单的 for 循环吗？首先，定义常量和函数。然后，运行 循环 n 次（示例中为 10'000），每次为变量绘制新的随机值并计算函数结果。最后，将所有结果附加到 results_dist，然后绘制它。
import numpy as np
import matplotlib.pyplot as plt
gym = [30, 30, 35, 35, 35, 35, 35, 35, 40, 40, 40, 45, 45]
brush_my_teeth = 2
size = 1000
def total_time(variables):
return np.sum(variables)
results_dist = []
for i in range(10000):
shower = np.random.triangular(left=5, mode=9, right=10, size)
argument = np.random.choice([0, 45], size, p=[0.9, 0.1])
dinner = np.random.normal(mu=15, sigma=5/3, size)
work = np.random.normal(mu=45, sigma=15/3, size)
variables = gym, shower, dinner, argument, work, brush_my_teeth
results_dist.append(total_time(variables))
plt.figure()
plt.hist(results_dist)
plt.show()
现有答案的想法是正确的，但我怀疑您想像 nicogen 那样对 size 中的所有值求和。
我假设您选择了一个相对较大的 size 来展示直方图中的形状，而不是您想要从每个类别中总结一个值。例如，我们要计算每个实例的一个实例的总和 activity，而不是 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]
你问的是获得"probability distribution"。像我们上面所做的那样获取直方图并不完全适合你。它为您提供了视觉表示，但不是分布函数。为了得到这个函数，我们需要使用核密度估计。 scikit-learn 有一个固定的 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
对于这类事情，我建议查看 Halton sequences and similar quasi-random low-discrepancy sequences. The ghalton 包，可以轻松生成确定性但差异较小的序列：
import ghalton as gh
sequence = gh.Halton(n) # n is the number of dimensions you want
然后根据其他一些答案，您可以执行以下操作：
values = sequence.get(10000) # generate a bunch of draws of
for vals in values:
# vals will have a single sample of n quasi-random numbers
variables = # add whatever other stuff you need to your quasi-random values
results_dist.append(total_time(variables))
如果您查看一些关于准随机序列的研究论文，它们已被证明可以更好地融合 Monte Carlo 集成和采样等应用程序。基本上，您可以更均匀地覆盖搜索 space，同时在样本中保持类似随机的属性，这在大多数情况下会导致更快的收敛。
这基本上让你在 n 维度上均匀分布。如果你想在某些维度上有非均匀分布，你可以相应地转换你的均匀分布。我不确定这会对 Halton 序列的低差异 属性 产生什么影响，但这可能值得研究。
给定一个依赖于多个变量的函数，每个变量都有一定的概率分布，我如何进行 Monte Carlo 分析以获得函数的概率分布。理想情况下，随着参数数量或迭代次数的增加，我希望解决方案具有高性能。
例如，我为 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)
健身房
淋浴
晚餐
参数
工作
brush_my_teeth
您尝试过简单的 for 循环吗？首先，定义常量和函数。然后，运行 循环 n 次（示例中为 10'000），每次为变量绘制新的随机值并计算函数结果。最后，将所有结果附加到 results_dist，然后绘制它。
import numpy as np
import matplotlib.pyplot as plt
gym = [30, 30, 35, 35, 35, 35, 35, 35, 40, 40, 40, 45, 45]
brush_my_teeth = 2
size = 1000
def total_time(variables):
return np.sum(variables)
results_dist = []
for i in range(10000):
shower = np.random.triangular(left=5, mode=9, right=10, size)
argument = np.random.choice([0, 45], size, p=[0.9, 0.1])
dinner = np.random.normal(mu=15, sigma=5/3, size)
work = np.random.normal(mu=45, sigma=15/3, size)
variables = gym, shower, dinner, argument, work, brush_my_teeth
results_dist.append(total_time(variables))
plt.figure()
plt.hist(results_dist)
plt.show()
现有答案的想法是正确的，但我怀疑您想像 nicogen 那样对 size 中的所有值求和。
我假设您选择了一个相对较大的 size 来展示直方图中的形状，而不是您想要从每个类别中总结一个值。例如，我们要计算每个实例的一个实例的总和 activity，而不是 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]
你问的是获得"probability distribution"。像我们上面所做的那样获取直方图并不完全适合你。它为您提供了视觉表示，但不是分布函数。为了得到这个函数，我们需要使用核密度估计。 scikit-learn 有一个固定的 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
对于这类事情，我建议查看 Halton sequences and similar quasi-random low-discrepancy sequences. The ghalton 包，可以轻松生成确定性但差异较小的序列：
import ghalton as gh
sequence = gh.Halton(n) # n is the number of dimensions you want
然后根据其他一些答案，您可以执行以下操作：
values = sequence.get(10000) # generate a bunch of draws of
for vals in values:
# vals will have a single sample of n quasi-random numbers
variables = # add whatever other stuff you need to your quasi-random values
results_dist.append(total_time(variables))
如果您查看一些关于准随机序列的研究论文，它们已被证明可以更好地融合 Monte Carlo 集成和采样等应用程序。基本上，您可以更均匀地覆盖搜索 space，同时在样本中保持类似随机的属性，这在大多数情况下会导致更快的收敛。
这基本上让你在 n 维度上均匀分布。如果你想在某些维度上有非均匀分布，你可以相应地转换你的均匀分布。我不确定这会对 Halton 序列的低差异 属性 产生什么影响，但这可能值得研究。