如何對方程進行蒙特卡羅分析？

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) 

健身房

淋浴

晚餐

參數

工作

brush_my_teeth

Answer 1

現有的答案有正確的想法，但我懷疑你要總結所有的值在size作爲nicogen已經完成。

我假設你是選擇一個比較大的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學會有一個罐頭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 

Answer 2

你有一個簡單的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() 

Answer 3

對於這樣的事情，我建議尋找到Halton sequences和類似準隨機低差異序列。該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積分和採樣等應用。基本上，您可以更均勻地覆蓋搜索空間，同時在樣本中保留隨機樣的屬性，這在大多數情況下會導致更快的收斂。

這基本上讓你在n尺寸上的統一分配。如果你想在某些維度上有不均勻的分佈，你可以相應地改變你的均勻分佈。我不確定這會對Halton序列的低差異性產生什麼影響，但它可能是值得研究的。