MM algorithm算法实现

算法来自于Abdallah, T. , & Vulcano, G. . (2020). Demand estimation under the multinomial logit model from sales transaction data. Manufacturing & Service Operations Management.

在python环境下，仿真得到了一些trsaction data，尝试实现算法，但不知道是否正确，求助如何正确实现以得到β和λ的估计

import pandas as pd
import numpy as np
import random
random.seed(42)

## initialize transaction data
# input data
n = 20
T = 100
delta = 0.9
la = 50

def get_trans(n=20,T=100,delta=0.9,la=50):
    '''
    n: number of items
    T: number of time period
    delta: probability that any item is available in any given period
    la: arrival rates of a homogenous poisson process
    '''
    ## simulate exp(beta)~Unif[0.05,1]
    exp_beta = [0]+[random.uniform(0.05,1) for i in range(n)]

    s = np.sum(exp_beta)/(1+np.sum(exp_beta)) # market share
    s_tilde = s/(1-s)

    arrival = np.random.poisson(lam=la, size=T)

    trans = np.zeros((n,T))
    for t in range(T):
        for i in range(n):
            if random.random() > delta:
                trans[i,t] = 0
            else:
                trans[i,t] = arrival[t]*exp_beta[i+1]*s
    return trans,s_tilde

def MM(trans,s_tilde):
    ## get n,T
    n = trans.shape[0]
    T = trans.shape[1]
    
    ## initialize beta
    lst = [0]+[np.log(1/n) for i in range(n)]
    Beta = np.array(lst)
    
    dist=100 # stopping criteria: beta间距离<0.1
    num = 0 # 统计迭代次数
    while dist >= 0.1 and num <= 10000:
        ### 计算q_it(Beta)
        q_beta = np.zeros((n,T))
        for t in range(T):
            sum_beta = np.sum(np.exp(Beta[np.where(trans[:,t]!=0)]))
            for i in range(n):
                if trans[i,t] == 0:
                    q_beta[i,t] = 0
                else:
                    q_beta[i,t] = np.exp(Beta[i+1])/sum_beta
        ### 计算beta
        for i in range(n):
            K_i = np.sum(trans[i,:])
            m_t = trans.sum(0)
            q_it = q_beta[i,:]
            mq_it = np.sum(m_t*q_it)
            Beta[i+1] = Beta[i+1] + np.log(K_i/mq_it)

        s_tilde1 = np.sum(np.exp(Beta))
        Beta_new = np.zeros(n+1)
        Beta_new[1:] = Beta[1:] + np.log(s_tilde/s_tilde1)
        dist = np.sum(np.abs(Beta_new - Beta))
        Beta = Beta_new
        num += 1
    return Beta

trans = get_trans()[0]
s_tilde = get_trans()[1]

Beta = MM(trans,s_tilde)
Beta

你好，我是有问必答小助手，非常抱歉，本次您提出的有问必答问题，技术专家团超时未为您做出解答

本次提问扣除的有问必答次数，将会以问答VIP体验卡（1次有问必答机会、商城购买实体图书享受95折优惠）的形式为您补发到账户。

​​​​因为有问必答VIP体验卡有效期仅有1天，您在需要使用的时候【私信】联系我，我会为您补发。