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

如何解决Python Statsmodels 中多项 Logistic 的最大似然估计

我试图模拟 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,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],1]})

y_full = np.ones(shape=(N*2,1))
class_1 = np.where(((U[:,1])),'class_1','class_0')
class_2 = np.where((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_classes])


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

Results from Scipy's minimize

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

解决方法

我发现了这个问题。因为多项 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]
    
    ## 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_classes])


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

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

The correct parameters found using scipy