Python统计参数估计代码优化求助：提升大重复次数运行速度

针对你的参数估计模拟代码，以下是几个关键优化点，能大幅提升运行速度：

1. 替换数值优化为解析解（核心优化） #

你的代码用pygosolnp.solve做数值优化来估计指数分布的scale参数，但指数分布的极大似然估计（MLE）有解析解：scale参数的估计值就是样本均值。这一步能直接消除最耗时的数值优化步骤，速度提升最明显。

指数分布的对数似然函数推导：

对于样本$v_1,...,v_n$，指数分布$Exp(scale=\alpha)$的对数似然为：
$L(\alpha) = \sum_{i=1}^n \left( -\ln\alpha - \frac{v_i}{\alpha} \right)$
对$\alpha$求导并令导数为0，解得$\hat{\alpha} = \frac{1}{n}\sum_{i=1}^n v_i$（样本均值）

所以完全不需要调用solve，直接用v.mean()就能得到参数估计值。

2. 避免循环中频繁调用numpy.append #

numpy.append每次都会创建新数组并复制数据，循环中多次调用会导致大量内存开销。应该预先分配数组空间，直接赋值。

3. 并行化重复实验 #

每个重复实验（RE次）是完全独立的，可以用多进程并行计算，利用CPU多核资源进一步加速。

优化后的完整代码 #

import numpy as np
from scipy.stats import skew, kurtosis
from tabulate import tabulate
from time import time
from multiprocessing import Pool

def simulate_single_run(args):
    n, alpha = args
    v = np.random.exponential(alpha, size=n)
    return v.mean()

def simulation(n_list, re, alpha):
    # 预先初始化统计结果数组
    stats = {
        "Mean": ["Mean"],
        "Variance": ["Variance"],
        "Bias": ["Bias"],
        "EQM": ["EQM"],
        "Skewness": ["Skewness"],
        "Kurtose": ["Kurtose"]
    }

    for n in n_list:
        # 并行计算所有重复实验的参数估计值
        with Pool() as pool:
            params = np.array(pool.map(simulate_single_run, [(n, alpha)] * re))
        
        # 计算统计量
        med = params.mean()
        varia = params.var()
        b = alpha - med
        eqma = b ** 2 + varia
        skewn = skew(params)
        kur = kurtosis(params)

        # 更新统计结果
        stats["Mean"].append(med)
        stats["Variance"].append(varia)
        stats["Bias"].append(b)
        stats["EQM"].append(eqma)
        stats["Skewness"].append(skewn)
        stats["Kurtose"].append(kur)

    # 整理输出表格
    data = list(stats.values())
    print(tabulate(data, headers=["Statistics", "n = 30", "n = 50", "n = 100", "n = 200", "n = 300"]))

if __name__ == "__main__":
    inicial_time = time()
    N = [30, 50, 100, 200, 300]
    RE = 100
    Alpha = 1/5
    simulation(N, RE, Alpha)
    print(f'运行时间: {time()-inicial_time} seconds')

额外优化建议 #

• 如果RE非常大（比如10000次），可以考虑用numpy向量化生成所有样本，一次性生成(re, n)的数组，然后直接按行求均值，避免多进程的开销：

  params = = np.random.exponential(Alpha, size=(RE, n)).mean(axis=1)
  

  这种方式在样本量不大时，速度可能比多进程更快。
• 减少不必要的导入，比如原代码中导入的log、expon.pdf在优化后都不需要了，可以删除。

内容的提问来源于stack exchange，提问作者TheJoing