如何在scipy.optimize.minimize中获取logit估计的准确标准误？

Great question—L-BFGS-B's Hessian approximation is indeed a known weak spot when you need reliable standard errors, since it's a limited-memory method that only tracks a small number of past updates instead of building the full Hessian. Here are several robust approaches to get accurate standard errors after your optimization:

1. Compute the Exact Hessian Directly #

Instead of relying on L-BFGS's approximate inverse Hessian, calculate the true Hessian (or the observed information matrix, which is the negative Hessian for MLE problems like logit) at the optimal parameter values. You can do this either analytically (if you can derive the Hessian formula) or numerically using a reliable differentiation library.

Example with numerical differentiation (using numdifftools):

import numdifftools as nd
import numpy as np

# Calculate the exact Hessian at the optimal parameters
hessian_exact = nd.Hessian(logit)(opt1.x)
# For logit MLE, standard errors are sqrt(diag(inverse(-Hessian)))
se_exact = np.sqrt(np.diag(np.linalg.inv(-hessian_exact)))

Note: Numerical differentiation can have precision issues for complex functions, but numdifftools handles this better than manual finite differences.

2. Switch to Full BFGS (If Memory Allows) #

If your parameter space isn't too large, use the full BFGS method instead of L-BFGS-B. BFGS maintains a complete approximation of the Hessian throughout optimization, which is far more reliable than L-BFGS's limited-memory version.

Example with SciPy's BFGS:

from scipy.optimize import minimize

opt_bfgs = minimize(logit, args=(df), x0=x_start, method='BFGS')
# BFGS's hess_inv is a full, more accurate approximation
hessinv_bfgs = opt_bfgs.hess_inv.todense()
se_bfgs = np.sqrt(np.diag(hessinv_bfgs))

If you need constrained optimization (the main reason to use L-BFGS-B), you can:

• Use a constrained method like SLSQP that supports better Hessian tracking
• Run L-BFGS-B to get the optimal parameters, then run a few iterations of unconstrained BFGS near that solution to refine the Hessian approximation.

3. Bootstrap for Non-Parametric Standard Errors #

Bootstrap is a robust, non-parametric approach that avoids Hessian calculations entirely. It works by resampling your dataset with replacement, refitting the model to each resample, and then computing the standard deviation of the parameter estimates across all resamples.

Example implementation:

def bootstrap_se(logit, df, x_start, n_iter=1000):
    n_samples = len(df)
    param_estimates = []
    
    for _ in range(n_iter):
        # Resample the dataset with replacement
        boot_df = df.sample(n_samples, replace=True)
        # Fit the model to the bootstrap sample
        opt = minimize(logit, args=(boot_df), x0=x_start, method='L-BFGS-B')
        param_estimates.append(opt.x)
    
    # Convert to array and compute standard deviations
    param_arr = np.array(param_estimates)
    se_bootstrap = np.std(param_arr, axis=0)
    return se_bootstrap

# Calculate bootstrap standard errors
se_boot = bootstrap_se(logit, df, x_start, n_iter=1000)

Tradeoff: Bootstrap is computationally expensive (especially with large datasets or many iterations) but produces highly reliable standard errors, even when Hessian-based methods fail.

4. Use the Score (Fisher) Information Matrix #

For maximum likelihood estimation (like your logit model), you can estimate the Fisher information using the outer product of score vectors instead of the Hessian. The score vector for each sample is the gradient of the log-likelihood for that sample, and the information matrix is the average outer product of these vectors.

Example for logit:

def logit_score(params, df):
    X = df.drop('y', axis=1).values
    y = df['y'].values
    y_hat = 1 / (1 + np.exp(-X @ params))
    # Score vector for each sample: (y_i - y_hat_i) * X_i
    scores = (y - y_hat)[:, np.newaxis] * X
    return scores

# Calculate score vectors at optimal parameters
scores = logit_score(opt1.x, df)
# Compute the observed score information matrix
info_matrix = np.dot(scores.T, scores) / len(df)
# Standard errors are sqrt(diag(inverse(info_matrix)))
se_score = np.sqrt(np.diag(np.linalg.inv(info_matrix)))

This method avoids calculating the Hessian entirely and is often more stable than numerical Hessian computation for MLE problems.

Each method has its pros and cons: exact Hessian is direct but may have precision issues; BFGS is accurate but memory-intensive; bootstrap is robust but slow; score information is efficient for MLE. Choose based on your dataset size, parameter count, and computational resources.

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