Python Statsmodels 中多项 Logistic 的最大似然估计答案

【问题标题】：Maximum Likelihood Estimate for Multinomial Logistic in Python's StatsmodelsPython Statsmodels 中多项 Logistic 的最大似然估计
【发布时间】：2021-05-27 06:57:03
【问题描述】：

我试图模拟 MNLogit 的 python statsmodels 实现只是为了更好地理解它，并且我可以重新创建与一些假数据上报告的分数相匹配的对数似然函数。

np.random.seed(100)
N = 1000
mu = [0,0]
rho = 0.1
cov = [[1, rho], [rho, 1]]

# u is N*2
u = np.random.multivariate_normal(mu, cov, 1000)
x1 = np.random.uniform(0, 1, size=(N,2)) #np.random.rand(N,2)
x2 = np.random.uniform(0, 1, size=(N,2)) #np.random.rand(N,2)

U = -1 + -3*x1 + 4*x2 + u

y = np.zeros(shape=(N, 2))
y[:,0] = ((U[:,0] > 0) & (U[:,0] > U[:,1]))
y[:,1] = (U[:,1] > 0 & (U[:,1] > U[:,0]))

W1 = pd.DataFrame({'x1':x1[:,0], 'x2':x2[:,0]})
W2 = pd.DataFrame({'x1':x1[:,1], 'x2':x2[:,1]})

y_full = np.ones(shape=(N*2,1))
class_1 = np.where(((U[:,0] > 0) & (U[:,0] > U[:,1])), 'class_1', 'class_0')
class_2 = np.where((U[:,1] > 0 & (U[:,1] > U[:,0])), 'class_2', 'class_0')
y_full = np.append(class_1, class_2)
W_full = sm.add_constant(W1.append(W2)).reset_index(drop=True)
mn_logit = sm.MNLogit(y_full, W_full)
mn_logit_res = mn_logit.fit()
mn_logit_res.summary()

这让我合理地拟合了我的假数据中的原始参数

MNLogit Results from statsmodels

如果我写出负似然函数并使用 scipy 最小化函数，我可以恢复相同的对数似然 (1260.8)，但参数估计值不同。

def cdf(W, beta):
    Wb = np.dot(W, beta)
    eXB = np.exp(Wb)
    eXB = eXB /eXB.sum(1)[:, None]
    return eXB

def take_log(probs):
    epsilon = 1e-20 
    return np.log(probs)

def calc_ll(logged, d):
    ll = d * logged
    return ll

def ll_mn_logistic(params, *args):
    y, W, n_params, n_classes = args[0], args[1], args[2][0], args[3]
    beta = [params[i] for i in range(0, len(params))]
    beta = np.array(beta).reshape(n_params, -1, order='F')
    
    ## onehot_encode
    d = pd.get_dummies(y, prefix='Flag').to_numpy()
    
    probs = cdf(W, beta)
    logged = take_log(probs)
    ll = calc_ll(logged, d)
    
    return -np.sum(ll)

n_params = 3, 
n_classes = 3
z = np.random.rand(3,3).flatten()
#probs = ll_mn_logistic(list(z), *[y_full, W_full, n_params, n_classes])


res = minimize(ll_mn_logistic, x0 =z, args =(y_full, W_full, n_params, n_classes), 
             options={'disp': True, 'maxiter':1000})
res

Results from Scipy's minimize

我怀疑这些差异是由于 statsmodels 中应用的优化方法造成的，但我已经尝试了 scipy 的一堆结果，但没有一个结果与使用 hessian 和 jacobian 的 statsmodels 返回的结果相近。谁能解释如何使用 scipy 的最小化来更好地近似 statsmodel 的方法？

【问题讨论】：

好问题！我建议您在问题中包含您正在谈论的具体价值观。我知道你包括截图。输出看起来不同，我不能很快明白你在说什么。现在，我会去看看我是否可以帮助找到答案或在哪里找到它。
看这里github.com/statsmodels/statsmodels/blob/…
你的 statsmodels 版本有 6 个参数，但你的有 9 个参数
两个版本的负对数似然优化值相同
感谢文森特，我已经更新了问题以包括在两者中看起来相同的负对数似然。 Josef 是的，我认为这是意料之中的。通常，通过估计 n-1 个回归并推断第三个回归来估计 n 个类别的 MLE。不完全确定为什么 statsmodel 摘要对象忽略报告第三类，但我输入了相同的数据，例如y_full, W_full 代表 3 个类。

标签： python scipy statsmodels scipy-optimize mle

【解决方案1】：

我发现了这个问题。因为多项式 MLE 算法适合除一个类回归之外的所有类回归，所以当我为三个类中的每一个传递参数时，我需要将 scipy 优化的 3 个参数集之一“归零”。

def cdf(W, beta):
    Wb = np.dot(W, beta)
    eXB = np.exp(Wb)
    eXB = eXB /eXB.sum(1)[:, None]
    return eXB

def take_log(probs):
    epsilon = 1e-20 
    return np.log(probs)

def calc_ll(logged, d):
    ll = d * logged
    return ll

def ll_mn_logistic(params, *args):
    y, W, n_params, n_classes = args[0], args[1], args[2][0], args[3]
    beta = [params[i] for i in range(0, len(params))]
    beta = np.array(beta).reshape(n_params, -1, order='F')

    #This step ensures that we're only fitting two parameter sets.
    beta[:,0] = [0, 0, 0]
    
    ## onehot_encode
    d = pd.get_dummies(y, prefix='Flag').to_numpy()
    
    probs = cdf(W, beta)
    logged = take_log(probs)
    ll = calc_ll(logged, d)
    
    return -np.sum(ll)

n_params = 3, 
n_classes = 3
z = np.random.rand(3,3).flatten()
#probs = ll_mn_logistic(list(z), *[y_full, W_full, n_params, n_classes])


res = minimize(ll_mn_logistic, x0 =z, args =(y_full, W_full, n_params, n_classes), 
             options={'disp': True, 'maxiter':1000})
res

这为我们提供了其他两个类的预期参数值。

The correct parameters found using scipy

【讨论】：