Python中的統計模擬

---

統計模擬是利用計算機方法從機率分佈中生成隨機樣本的任務，以便我們可以對錶現出隨機行為的複雜系統進行建模和分析。在本文中，我們將瞭解如何在Python中使用這個強大的工具進行預測、生成見解以及評估統計算法的效能。

存在不同型別的統計模擬，如下所示：

• 蒙特卡洛模擬 - 從機率分佈生成隨機樣本以估計函式的期望值。

• Bootstrap - 常用於估計估計量抽樣分佈的重抽樣技術。

• 馬爾可夫鏈蒙特卡洛 (MCMC) - 用於估計複雜機率分佈引數的一類演算法。

• 隨機過程模擬 - 模擬隨時間變化的隨機行為。股票價格或天氣模式預測是一些例子。

統計模擬廣泛應用於金融、工程、物理、生物和社會科學等領域，用於對複雜系統進行建模和分析、進行預測和生成見解。使用基於計算機的方法可以有效且準確地模擬大量場景，這些場景可用於量化不確定性、估計風險並做出資料驅動的決策。

蒙特卡洛模擬

蒙特卡洛模擬是一種統計模擬，它涉及從機率分佈生成隨機樣本以估計函式的期望值，或對錶現出隨機行為的複雜系統進行建模和分析。“蒙特卡洛”這個名稱源於位於摩納哥的蒙特卡洛賭場，在那裡，隨機性是機會遊戲中，更準確地說是在賭博中的關鍵方面。

蒙特卡洛模擬的精度取決於生成的隨機樣本數量以及被分析的模型或系統的質量。如果樣本數量足夠多且模型良好，蒙特卡洛模擬可以提供寶貴的見解，並幫助決策者做出明智的選擇。

示例

import numpy as np

# Define the function to be evaluated
def function(x):
   return x**2

# Generate random samples from a uniform distribution between 0 and 1
samples = np.random.uniform(0, 1, size=10000)

# Evaluate the function at each sample
values = function(samples)

# Compute the average of function values
mean_value = np.mean(values)

print("Mean value of the function: ", mean_value)

輸出

Mean value of the function:  0.3326046914715845

上面的程式碼演示了一個非常簡單的示例，說明我們如何使用蒙特卡洛模擬來估計函式的平均值。透過生成大量隨機樣本，然後在每個樣本處評估函式，我們可以輕鬆獲得平均值的近似值。估計的精度將在很大程度上取決於生成的樣本數量以及我們將要評估的函式的複雜性。

Bootstrap

Bootstrap 是一種統計方法，用於透過對資料進行有放回的重抽樣來估計估計量的抽樣分佈。這是一種強大的技術，可用於估計估計量的可變性並構建置信區間。該方法最初由 Bradley Efron 先生於 1979 年提出。

當樣本量較小或總體分佈未知或複雜且中心極限定理等傳統方法不適用時，Bootstrap 特別有用。

Bootstrap 的基本步驟是：

• 收集資料樣本。

• 從原始樣本中抽取大量 Bootstrap 樣本（有放回）。

• 計算每個 Bootstrap 樣本感興趣的估計量（例如，均值、標準差等）。

• 使用從 Bootstrap 樣本計算出的估計量分佈來推斷總體。

Bootstrap 可應用於各種估計量，包括均值、標準差、迴歸係數等等。它還可以用於估計檢驗統計量的分佈，例如 t 統計量，可用於檢驗關於總體引數的假設。

以下是如何使用 Bootstrap 方法估計資料樣本的標準差並構建 95% 置信區間的示例。

示例

import numpy as np

# original sample
data = [1, 2, 3, 4, 5]

# number of bootstrap samples
n_samples = 1000

# array to store the bootstrap samples standard deviation
std_devs = np.empty(n_samples)

# generate bootstrap samples
for i in range(n_samples):
   sample = np.random.choice(data, size=len(data), replace=True)
   std_devs[i] = np.std(sample)

# calculating the lower as well as upper bound of the confidence interval
alpha = 0.05
lower = np.percentile(std_devs, alpha/2*100)
upper = np.percentile(std_devs, (1-alpha/2)*100)

print(f'Confidence interval: [{lower}, {upper}]')

輸出

Confidence interval: [0.4898979485566356, 1.7435595774162693]

在此示例中，我們從原始樣本中抽取 1000 個 Bootstrap 樣本並計算每個樣本的標準差。然後，我們使用標準差的分佈來計算 95% 置信區間的下限和上限，方法是使用 np.percentile() 函式。

馬爾可夫鏈蒙特卡洛 (MCMC)

馬爾可夫鏈蒙特卡洛 (MCMC) 是一類用於估計複雜機率分佈引數的演算法。這些演算法構建一個馬爾可夫鏈，該鏈具有所需的機率分佈作為其平衡分佈。透過模擬該鏈，可以從該分佈中生成樣本，然後將其用於統計推斷。

MCMC 的基本思想是構建一個馬爾可夫鏈，該鏈具有目標分佈作為其平穩分佈。然後執行該鏈許多步驟，並定期收集樣本。然後使用這些樣本估計目標分佈的引數。

以下是如何使用 Metropolis-Hastings 演算法（MCMC 的子部分）從正態分佈中抽樣的示例。

示例

import numpy as np

# target distribution
mean = 0
std = 1
target = lambda x: 1/(std * np.sqrt(2 * np.pi)) * np.exp(-(x - mean)**2 / (2 * std**2))

# proposal distribution
proposal_mean = 0
proposal_std = 2
proposal = lambda x: 1/(proposal_std * np.sqrt(2 * np.pi)) * np.exp(-(x - proposal_mean)**2 / (2 * proposal_std**2))

# initial state
x = 0

# number of samples
n_samples = 10000

# array to store the samples
samples = [x]

# Metropolis-Hastings algorithm
for i in range(n_samples):
   x_new = np.random.normal(x, proposal_std)
   acceptance_prob = min(1, target(x_new) / target(x))
   if np.random.rand() < acceptance_prob:
      x = x_new
   samples.append(x)

# plot the samples
import matplotlib.pyplot as plt
plt.hist(samples, bins=50, density=True)
plt.show()

輸出

此指令碼使用 Metropolis-Hastings 演算法從均值為 0 且標準差為 1 的正態分佈生成一系列樣本，使用均值為 0 且標準差為 2 的正態建議分佈。該演算法從初始狀態開始（在本例中為 0），然後透過從建議分佈中提出新值來生成新狀態。根據由接受機率確定的機率接受新值，該機率是根據目標分佈和建議分佈的比率計算的。

隨機過程模擬

隨機過程模擬是一類統計模擬，涉及模擬隨時間變化的隨機行為。它們用於對錶現出隨機性的複雜系統進行建模和分析，例如股票價格、天氣模式和生物種群。

隨機過程是一種數學模型，它描述了一個其行為受隨機性影響的系統。

示例

import numpy as np
import matplotlib.pyplot as plt

# Parameters
p = 0.5  # probability of coin being a head
T = 10  # defining the number of time steps

# Initial state
x = 1

# Random number generator
np.random.seed(0)

# Array to store the states
states = [x]

# Simple stochastic process simulation
for t in range(T):
    x = 1 if np.random.rand() < p else 0
    states.append(x)

print(states)

# Plot the states
plt.step(range(T+1), states, where='post')
plt.xlabel('Time')
plt.ylabel('State')
plt.ylim(-0.5,1.5)
plt.show()

輸出

[1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1]

其他基本示例

以下是如何模擬擲骰子的示例。

示例

import numpy as np

# Generate random samples from a uniform distribution between 1 and 6
samples = np.random.randint(1, 7, size=10000)

# Compute the average and the standard deviation of the samples
sample_mean = np.mean(samples)
sample_std = np.std(samples)

print("Sample mean: ", sample_mean)
print("Sample standard deviation: ", sample_std)

輸出

Sample mean:  3.4946
Sample standard deviation:  1.7094358250604205

結論

統計模擬是一個強大的工具，可以幫助理解和分析複雜的系統和過程。Python 為我們提供了一系列工具和庫，使統計模擬的實現和執行變得容易。

統計模擬的主要優點是它不僅能夠對複雜的系統和過程進行建模，而且能夠研究不同引數和變數對結果的影響。它的應用幾乎對所有領域都有幫助。

NumPy、SciPy 和 Pandas 是一些可用於生成隨機樣本、評估函式和執行統計分析的庫示例。

總而言之，統計模擬對於從事資料科學或相關領域工作的任何人來說都是必不可少的工具，而 Python 為實現這些模擬提供了一個靈活且強大的平臺。透過利用統計模擬的強大功能，我們可以獲得對複雜系統和過程的新見解，並根據我們的分析做出更明智的決策。