如何用Maximum Likelihood Estimator估计幂律功率谱的指数？

核心需求 #

已知离散功率谱满足幂律分布 f(x) = A·x^b，需要用最大似然估计（MLE）拟合指数b，规避双对数最小二乘方法的固有偏差。

实现原理 #

功率谱的观测值通常服从伽马分布，对应的负对数似然函数为：

L(b, A) = Σ(ln(A·x^b) + f(x)/(A·x^b))
我们需要最小化这个负对数似然来得到最优的b和A。注意参数A必须为正，优化时需添加边界约束避免无效值。

完整代码实现 #

from astroML.time_series import generate_power_law
import numpy as np
import scipy
from scipy.optimize import minimize

# 生成模拟数据
N = 72  # 数据点数量
dt = 5 * 60  # 时间分辨率
true_exponent = -2  # 真实幂律指数
rand_seed = 1
time_max = 6 * 6 * 60  # 观测时长

# 生成符合幂律的时间序列
y = generate_power_law(N, dt, -true_exponent, random_state=rand_seed) 

# 计算功率谱密度（PSD）
sample_freq = scipy.fft.rfftfreq(N, d=dt)[1:]  # 采样频率（去除直流分量）
sig_fft = scipy.fft.rfft(y, N)[1:]  # FFT结果（去除直流分量）
psd = (np.abs(sig_fft)**2) * 2 * dt / N  # 功率谱密度

# 定义幂律模型
def powerlaw_model(x, b, A):
    return A * (x ** b)

# 定义负对数似然函数
def neg_log_likelihood(params, x, y):
    b, log_A = params
    A = np.exp(log_A)  # 对数变换确保A始终为正
    model = powerlaw_model(x, b, A)
    # 避免模型值为0/负数导致计算错误
    if np.any(model <= 0):
        return 1e30
    return np.sum(np.log(model) + y / model)

# 初始参数猜测
initial_guess = [-1, 0]  # b初始设为-1，log(A)初始设为0（对应A=1）

# 优化：选用支持边界约束的L-BFGS-B方法
result = minimize(neg_log_likelihood, initial_guess, args=(sample_freq, psd),
                  method='L-BFGS-B')

# 提取并打印结果
estimated_b = result.x[0]
estimated_A = np.exp(result.x[1])

# 计算最小二乘结果用于对比
log_freq = np.log(sample_freq)
log_psd = np.log(psd)
ls_slope, _ = np.polyfit(log_freq, log_psd, 1)

print(f"真实指数: {true_exponent}")
print(f"MLE估计指数: {estimated_b:.4f}")
print(f"MLE估计系数A: {estimated_A:.4f}")
print(f"最小二乘拟合指数: {ls_slope:.4f}")

关键说明 #

• 参数约束：通过对A做对数变换，确保优化过程中A始终为正，避免出现无效的负系数。
• 优化方法选择：使用L-BFGS-B而非Nelder-Mead，前者支持边界约束，收敛性和稳定性更优。
• 偏差消除：MLE估计的指数会比双对数最小二乘结果更接近真实值，彻底解决最小二乘的固有偏差问题。

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