二项分布对数MLE模拟遇RuntimeWarning：log除零错误排查

问题背景 #

我在模拟样本量增加时伯努利分布（二项分布n=1的情况）对数极大似然估计（log MLE）的变化，真实参数设为0.6。已经排除了θ=0的情况，但代码第19行还是出现RuntimeWarning: divide by zero encountered in log错误，错误行是cal_llh = np.log(theta**(x) * (1-theta)**(1-x))。

原始代码 #

############################################################
## Step 1#
############################################################
# function to calculate likelihood for bernoulli
def logLikelihood(theta, x):
    # cal the log likelihood of each observation in the samples collected
    cal_llh = np.log(theta**(x) * (1-theta)**(1-x))
    tlld = np.prod(cal_llh)# cal the total likleihood
    return tlld

# function to calculate
def mle_Binom(X_samples, thetas):
    loglikelihood_single_theta = [logLikelihood(theta=t, x=X_samples) for t in thetas]
    # mle_val=thetas[np.argmax(likelihood_single_theta)] #get the maximum likelihood estimate
    return np.array(loglikelihood_single_theta)

# test the functions
true_params_Bern = 0.6


############################################################
## Step 2#
############################################################
# how does the likelihood plot changes as sample size changes 
Bern_Nsamples = np.linspace(start=100, stop=1000, num=100, dtype=int)
response_Bernoulli = np.random.binomial(n=1, p=0.6, size=100)
possible_thetas = np.linspace(start=0.001, stop=1, num=100)
result_theta = np.ma.array(possible_thetas.copy())

Bern_Nsamples = np.linspace(start=100, stop=1000, num=100, dtype=int)
beta_for_mle_holder = []

def Bernoulli_optim_nSamples(Bern_stepSize, rand_sets_thetas): 
    for n in Bern_stepSize:
        response_Bernoulli = np.random.binomial(n=1, p=0.6, size=n)
        mle_out_Binom = mle_Binom(X_samples=response_Bernoulli, thetas=rand_sets_thetas) #cal lld of specific theta
        max_theta_Binom = rand_sets_thetas[np.argmax(mle_out_Binom)] #which theta gave us max lld
        beta_for_mle_holder.append(max_theta_Binom)
    fig, ax = plt.subplots()
    ax.plot(Bern_stepSize, beta_for_mle_holder)
    ax.set_title('Bernoulli dist nSamples vrs MLE')
    ax.hlines(y=0.6, xmin=min(Bern_stepSize), xmax=max(Bern_stepSize), linestyles="dashed", color="red", label="MLE")
    plt.xlabel("nSamples")
    plt.ylabel("MLE")
    plt.show()


Bernoulli_optim_nSamples(Bern_stepSize=Bern_Nsamples, rand_sets_thetas=result_theta)

错误信息 #

fods23/simulations/scripts/Binomial_MLE_simulations.py:19: RuntimeWarning: divide by zero encountered in log
  cal_llh = np.log(theta**(x) * (1-theta)**(1-x))

错误原因 #

1. θ取值包含1导致log(0)：你的possible_thetas包含了θ=1的情况，当θ=1时，1-theta=0。如果样本中存在x=0的观测值，(1-theta)**(1-x)会变成0^1=0，最终theta**x*(1-theta)**(1-x)的结果为0，取log时就会触发除零警告（因为log(0)在数值计算中等价于负无穷，会被判定为除零错误）。
2. 对数似然计算逻辑错误：对数似然的总和应该是对每个观测的对数似然求和，而不是取乘积。原始似然是各观测似然的乘积，取log后转化为求和，你当前用np.prod(cal_llh)完全搞反了逻辑，会导致数值结果错误甚至溢出。

修复方案 #

1. 调整θ的取值范围 #

将θ的上限从1改为0.999，避免出现1-theta=0的情况：

possible_thetas = np.linspace(start=0.001, stop=0.999, num=100)

2. 修正对数似然计算逻辑 #

利用对数性质拆分计算，避免先算乘积再取log（减少数值下溢风险），同时将乘积改为求和：

def logLikelihood(theta, x):
    # 利用对数性质拆分计算，避免直接计算乘积
    cal_llh = x * np.log(theta) + (1 - x) * np.log(1 - theta)
    # 对数似然总和是求和，不是乘积
    tlld = np.sum(cal_llh)
    return tlld

3. 优化代码效率（可选） #

用numpy向量化操作代替列表推导，提升计算速度：

def mle_Binom(X_samples, thetas):
    # 向量化计算所有theta的对数似然
    x = X_samples[:, np.newaxis]
    loglikelihoods = x * np.log(thetas) + (1 - x) * np.log(1 - thetas)
    return loglikelihoods.sum(axis=0)

完整修正后代码 #

import numpy as np
import matplotlib.pyplot as plt

############################################################
## Step 1
############################################################
# 修正后的伯努利分布对数似然计算函数
def logLikelihood(theta, x):
    # 利用对数性质拆分，避免数值下溢
    cal_llh = x * np.log(theta) + (1 - x) * np.log(1 - theta)
    # 对数似然总和为求和
    total_log_likelihood = np.sum(cal_llh)
    return total_log_likelihood

# 向量化优化的MLE计算函数
def mle_Binom(X_samples, thetas):
    # 将样本转为列向量，与theta的行向量广播计算
    x_col = X_samples.reshape(-1, 1)
    loglikelihoods = x_col * np.log(thetas) + (1 - x_col) * np.log(1 - thetas)
    return loglikelihoods.sum(axis=0)

# 真实参数
true_params_Bern = 0.6


############################################################
## Step 2
############################################################
# 样本量序列
Bern_Nsamples = np.linspace(start=100, stop=1000, num=100, dtype=int)
# 排除0和1的theta取值范围
possible_thetas = np.linspace(start=0.001, stop=0.999, num=100)
result_theta = possible_thetas.copy()

beta_for_mle_holder = []

def Bernoulli_optim_nSamples(Bern_stepSize, rand_sets_thetas): 
    for n in Bern_stepSize:
        # 生成对应样本量的伯努利样本
        response_Bernoulli = np.random.binomial(n=1, p=0.6, size=n)
        # 计算所有theta的对数似然
        mle_out_Binom = mle_Binom(X_samples=response_Bernoulli, thetas=rand_sets_thetas)
        # 找到对数似然最大的theta
        max_theta_Binom = rand_sets_thetas[np.argmax(mle_out_Binom)]
        beta_for_mle_holder.append(max_theta_Binom)
    
    # 绘图
    fig, ax = plt.subplots()
    ax.plot(Bern_stepSize, beta_for_mle_holder, label="Estimated MLE")
    ax.set_title('Bernoulli Distribution: Sample Size vs MLE')
    ax.hlines(y=0.6, xmin=min(Bern_stepSize), xmax=max(Bern_stepSize), linestyles="dashed", color="red", label="True Parameter")
    plt.xlabel("Number of Samples")
    plt.ylabel("MLE of p")
    plt.legend()
    plt.show()


# 运行模拟
Bernoulli_optim_nSamples(Bern_stepSize=Bern_Nsamples, rand_sets_thetas=result_theta)

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