我怎样才能运行 5000 个样本的代码而不会永远花费时间

Question

我正在使用 5000 个样本（每个样本有 100 个数据点）来估计 beta 分布的最大似然估计量 (MLE)。然而，我目前的实现需要很长时间才能运行。

这是我正在使用的代码：

def U_score(x, theta):
    a = theta[0]
    b = theta[1]
    n = len(x)
    epsilon = 1e-10 # to avoid log(0)
    d_a = -sp.digamma(a) + sp.digamma(a + b) + sum(np.log(x + epsilon)) / (n)
    d_b = -sp.digamma(b) + sp.digamma(a + b) + sum(np.log(1 - x + epsilon)) / (n)
    return np.array((d_a, d_b))

def Hessiana(x, theta):
    a = theta[0]
    b = theta[1]
    h11 = sp.polygamma(1, a + b) - sp.polygamma(1, a)
    h12 = sp.polygamma(1, a + b)
    h21 = sp.polygamma(1, a + b) 
    h22 = sp.polygamma(1, a + b) - sp.polygamma(1, b)
    return np.array([[h11, h12], [h21, h22]])

def H_inv(x, theta):
    H = Hessiana(x, theta)
    ridge = 1e-6  # Small constant
    H_ridge = H + ridge * np.eye(H.shape[0])
    return np.linalg.inv(H_ridge)

def max_likelihood(x, theta, tol):
    iter = 0
    while iter < 1000:
        theta_new = theta - H_inv(x, theta) @ U_score(x, theta)
        if np.linalg.norm(theta_new - theta) <= tol:
            return theta_new, iter + 1
        theta = theta_new
        iter += 1
    print('Não convergiu!')
    return theta_new, iter

def mse(arr, theta_real):
    a = theta_real[0]
    b = theta_real[1]
    eqm_a = np.mean((arr[:,0] - a)**2)
    eqm_b = np.mean((arr[:,1] - b)**2)
    return mse_a, mse_b

arr = np.array([])
n_amostras = 5000
theta = np.array([1,1])
for i in range(0, n_amostras):
    x = st.beta.rvs(2, 5, size = 100)
    emv, iteracoes = max_likelihood(x, theta, 1e-6)
    arr = np.append(arr, emv)
arr = arr.reshape(-1, 2)
print(arr)

我需要计算的最后一部分是问题：

arr = np.array([])
n_amostras = 5000
theta = np.array([1,1])
for i in range(0, n_amostras):
    x = st.beta.rvs(2, 5, size = 100)
    emv, iteracoes = max_likelihood(x, theta, 1e-6)
    arr = np.append(arr, emv)
arr = arr.reshape(-1, 2)
print(arr)

如果您能提供有关如何优化此过程并更有效地计算均方误差 (MSE) 的建议，我将不胜感激。谢谢！

Answer 1

1. 为了更有效地计算均方误差（MSE），我建议您阅读有关堆栈溢出流的问题Numpy 中的均方误差？;

2. 不用python，你可以使用Pypy，那就是

   jit

   ，不需要对代码做任何改动；

3. 我尝试在不更改逻辑的情况下重写代码，并设法节省了大约 120,000 次函数调用（大约 3,000,000 次）。

代码下方：

from numpy import array, log, sum as sumnp, append
from numpy.linalg import norm, inv
from scipy.stats import beta
from scipy.special import digamma, polygamma

def U_score(x, theta):
    a, b = theta[0], theta[1]
    n = len(x)
    dig_ab = digamma(a + b)
    epsilon = 1e-10
    d_a = - digamma(a) + dig_ab + sumnp(log(x + epsilon)) / n
    d_b = - digamma(b) + dig_ab + sumnp(log(1 - x + epsilon)) / n
    return array((d_a, d_b))

def Hessiana(x, theta):
    a, b = theta[0], theta[1]
    pol_ab= polygamma(1, a + b)
    return array([[pol_ab - polygamma(1, a), pol_ab], [pol_ab, pol_ab - polygamma(1, b)]])

    def H_inv(x, theta):
    H = Hessiana(x, theta)
    return inv(H + 1e-6 * np.eye(H.shape[0]))

def max_likelihood(x, theta, tol):
    for cycle in range(1000):
        theta_new = theta - H_inv(x, theta) @ U_score(x, theta)
        if norm(theta_new - theta) <= tol:
            return theta_new, cycle + 1
        theta = theta_new
    raise Exception('Não convergiu!')


arr = array([])
n_amostras = 5000
theta = array([1, 1])
for i in range(n_amostras):
    x = beta.rvs(2, 5, size=100)
    emv, iterazioni = max_likelihood(x, theta, 1e-6)
    arr = append(arr, emv)
arr = arr.reshape(-1, 2)
print(arr)