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 为实现这些模拟提供了一个灵活且强大的平台。通过利用统计模拟的强大功能，我们可以获得对复杂系统和过程的新见解，并根据我们的分析做出更明智的决策。