First, let's load the libraries we shall need.
import numpy as np
import os
import pandas as pd
import string as str
import math
import sys
import time
import scipy.sparse as spr
from scipy import optimize, special
import gurobipy as grb
We will go back to the dataset of Greene and Hensher (1997). As a reminder, 210 individuals are surveyed about their choice of travel mode between Sydney, Canberra and Melbourne, and the various costs (time and money) associated with each alternative. Therefore there are 840 = 4 x 210 observations, which we can stack into travelmodedataset
a 3 dimensional array whose dimensions are mode,individual,dummy for choice+covariates.
Let's load the dataset and represent it conveniently in a similar fashion as in block 6:
thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/demand_travelmode/'
travelmode = pd.read_csv(thepath+'travelmodedata.csv')
travelmode['choice'] = np.where(travelmode['choice'] =='yes' , 1, 0)
nobs = travelmode.shape[0]
ncols = travelmode.shape[1]
nbchoices = 4
ninds = int(nobs/nbchoices)
muhat_i_y = travelmode['choice'].to_numpy().reshape(ninds,nbchoices).T
muhat_iy = muhat_i_y.flatten()
muhat_i_y = travelmode['choice'].to_numpy().reshape(ninds,4).T
muhat_iy = muhat_i_y.flatten()
s_y = travelmode.groupby(['mode']).mean()['choice'].to_frame().sort_index()
def two_d(X):
return np.reshape(X,(X.size, 1))
Start with assuming that there is no observable heterogeneity, so the only observation we have at hand are the aggregate market shares $s_y$. Hence the systematic utility will be the same for every agent. However, we wish to write a parametric model for it, namely assume a knwon parametric form for the dependence of $U_y$ with respect to various observed characteristics associated with $y$.
Assume then that the utilities are parameterized as follows: $U = \Phi \beta$ where $\beta\in\mathbb{R}^{p}$ is a parameter, and $\Phi$ is a $\left\vert \mathcal{Y}\right\vert \times p$ matrix.
The log-likelihood function is given by
\begin{align*} l\left( \beta\right) =N\sum_{y}\hat{s}_{y}\log\sigma_{y}\left(\Phi \beta\right) \end{align*}A common estimation method of $\beta$ is by maximum likelihood%
\begin{align*} \max_{\beta}l\left( \beta\right) . \end{align*}MLE is statistically efficient; the problem is that the problem is not guaranteed to be convex, so there may be computational difficulties (e.g. local optima).
In the logit case,
\begin{align*} l\left( \beta\right) =N\left\{ \hat{s}^{\intercal}\Phi\beta-\log\sum_{y}\exp\left( \Phi\beta\right) _{y}\right\} \end{align*}so that the max-likehood amounts to
\begin{align*} \max_{\beta}\left\{ \hat{s}^{\intercal} \Phi \beta-G\left( \Phi \beta\right) _{y}\right\} \end{align*}whose value is the Legendre-Fenchel transform of $\beta\rightarrow G\left( \Phi \beta\right)$ evaluated at $\Phi ^{^{\intercal}}\hat{s}$.
Note that the vector $\Phi^{^{\intercal}}\hat{s}$ is the vector of empirical moments, which is a sufficient statistics in the logit model.
As a result, in the logit case, the MLE is a convex optimization problem, and it is therefore both statistically efficient and computationally efficient.
The previous remark will inspire an alternative procedure based on the moments statistics $\Phi^{^{\intercal}}\hat{s}$.
The social welfare is given in general by $W\left( \beta\right) =G\left( \Phi\beta\right) $. One has $\partial_{\beta^{i}}W\left(\beta\right) =\sum_{y}\sigma_{y}\left( \Phi\beta\right) \Phi_{yi}$, that is
\begin{align*} \nabla W\left( \beta\right) = \Phi^{\intercal}\sigma\left( \Phi\beta\right) , \end{align*}which is the vector of predicted moments.
Therefore the program
\begin{align*} \max_{\beta}\left\{ \hat{s}^{\intercal}\Phi\beta-G\left( \Phi\beta\right) _{y}\right\} \end{align*}picks up the parameter $\beta$ which matches the empirical moments $X^{^{\intercal}}\hat{s}$ with the predicted ones $\nabla W\left(\beta\right) $. This procedure is not statistically efficient, but is computationally efficient becauses it arises from a convex optimization problem.
Back to the logit case. Recall we have
\begin{align*} l\left( \beta\right) =N\left\{ \hat{s}^{\intercal}\Phi\beta-\log\sum_{y} \exp\left( \Phi\beta\right) _{y}\right\} \end{align*}Assume that we restrict ourselves to $\beta^{\top}z>0$. Then we can write $\beta=\theta/T$ where $T=1/\beta^{\top}z$ and $\theta=\beta T$. Call $\Theta=\left\{ \theta\in\mathbb{R}^{p},\theta^{\top}z=1\right\} $, so that $\beta=\theta/T$ where $\theta\in\Theta$ and $T>0$. We have
\begin{align*} l\left( \theta,T\right) =\frac{N}{T}\left\{ \hat{s}^{\intercal} \Phi\theta-T\log\sum_{y}\exp\left( \frac{\left( \Phi\theta\right) _{y}}{T}\right) \right\} \end{align*}and we define the fixed temperature maximum likelihood estimator by
\begin{align*} \theta\left( T\right) =\arg\max_{\theta}l\left( \theta,T\right) \end{align*}Note that $\theta\left( T\right) =\arg\max_{\theta\in\Theta}Tl\left(\theta,T\right)$ where
\begin{align*} Tl\left( \theta,T\right) =N\left\{ \hat{s}^{\intercal}\Phi\theta-T\log\sum _{y}\exp\left( \frac{\left( \Phi\theta\right) _{y}}{T}\right) \right\} \end{align*}and we note that $Tl\left( \theta,T\right) \rightarrow N\left\{ \hat{s}^{\intercal}\Phi\theta-\max_{y\in\mathcal{Y}}\left\{ \left( \Phi\theta\right)_{y}\right\} \right\} $ as $T\rightarrow0$.
We have
\begin{align*} \frac{Tl\left( \theta,T\right) }{N}=\hat{s}^{\intercal}\Phi\theta-T\log\sum_{y}\exp\left( \frac{\left( \Phi\theta\right) _{y}}{T}\right) \end{align*}Let $\theta\left( 0\right) =\lim_{T\rightarrow0}\theta\left(T\right) $. Calling $m\left( \theta\right) =\max_{y\in\mathcal{Y}}\left\{\left( \Phi\theta\right) _{y}\right\} $, we have
\begin{align*} \theta\left( 0\right) \in\arg\max_{\theta}\left\{ \hat{s}^{\intercal}\Phi\theta-m\left( \theta\right) \right\}, \end{align*}or
\begin{align*} \theta\left( 0\right) \in\arg\min_{\theta}\left\{ m\left( \theta\right)-\hat{s}^{\intercal}\Phi\theta\right\}, \end{align*}Calling $m\left( \theta\right) =\max_{y\in\mathcal{Y}}\left\{ \left(\Phi\theta\right) _{y}\right\} $, one has
\begin{align*} \theta\left( T\right) \in\arg\max\left\{ \hat{s}^{\intercal}\Phi\theta-m\left( \theta\right) -T\log\sum_{y}\exp\left( \frac{\left(\Phi\theta\right) _{y}-m\left( \theta\right) }{T}\right) \right\} \end{align*}Note that
\begin{align*} \theta\left( 0\right) \in\arg\max\left\{ \hat{s}^{\intercal}\Phi\theta -m\left( \theta\right) \right\} . \end{align*}Define $R_{i}\left( \theta,y\right) =\left( \Phi\theta\right)_{y}-\left( \Phi\theta\right) _{y_{i}}$ the regret associated with observation $i$ with respect to $y$. This is equal to the difference between the payoff given by $y$ and the payoff obtained under observation $i$, denoting $y_{i}$ the action taken in observation $i$. The max-regret associated with observation $i$ is therefore
\begin{align*} \max_{y\in\mathcal{Y}}R_{i}\left( \theta,y\right) =\max_{y\in\mathcal{Y}}\left\{ \left( \Phi\theta\right)_{y}-\left( \Phi\theta\right)_{y_{i}}\right\} \end{align*}and the max-regret associated with the sample is $\frac{1}{N}\sum\max_{y\in\mathcal{Y}}\left\{ R_{i}\left( \theta,y\right) \right\} $, that is $\max_{y\in\mathcal{Y}}\left\{ \left( \Phi\theta\right) _{y}\right\} - \hat{s}^{\intercal}X\theta$.
The minimax regret estimator
\begin{align*} \hat{\theta}^{MMR}=\min_{\theta}\left\{ m\left( \theta\right) -\hat {s}^{\intercal}\Phi\theta\right\} \end{align*}which has a linear programming fomulation
\begin{align*} & \min_{m,\theta}m-\hat{s}^{\intercal}\Phi\theta\\ s.t.~ & m-\left( \Phi\theta\right) _{y}\geq\forall y\in\mathcal{Y} \end{align*}Note that the set of $\theta$ that enter the solution to the problem above is not unique, but is a convex set. Denoting $V$ the value of program, we can look for bounds of $\theta^{\intercal}d$ for a chosen direction $d$ by
\begin{align*} & \min_{\theta,m}/\max_{\theta,m} \theta^{\intercal}d\\ s.t.~ & m-\hat{s}^{\intercal}X\theta=V\\ & m\geq\left( \Phi\theta\right)_{y}, \quad \forall y\in\mathcal{Y}% \end{align*}See class notes
We now assume that we observe individual characteristics that are relevant for individual choices, that is $U_{iy}=\sum_k \Phi_{iyk} \beta_k$, or in matrix form $$U = \Phi \beta,$$ where $\beta\in\mathbb{R}^{p}$ is a parameter, and $\Phi$ is a $\left(\left\vert \mathcal{I}\left\vert\right\vert\mathcal{Y}\right\vert \right) \times p$ matrix.
See class notes.
Back to the dataset:
Phi_iy_k = np.column_stack((np.kron(np.identity(4)[0:4,1:4],np.repeat(1, ninds).reshape(ninds,1)), - travelmode['travel'].to_numpy(), - (travelmode['travel']*travelmode['income']).to_numpy(), - travelmode['gcost'].to_numpy()))
nbK = Phi_iy_k.shape[1]
phi_mean = Phi_iy_k.mean(axis = 0)
phi_stdev = Phi_iy_k.std(axis = 0, ddof = 1)
Phi_iy_k = ((Phi_iy_k - phi_mean).T/phi_stdev[:,None]).T
def log_likelihood(theta):
nbK = np.asarray(theta).shape[0]
Xtheta = Phi_iy_k.dot(theta)/sigma
Xthetamat_iy = Xtheta.reshape(nbchoices, ninds).T
max_i = np.amax(Xthetamat_iy, axis = 1)
expPhi_iy = np.exp((Xthetamat_iy.T -max_i).T)
d_i = np.sum(expPhi_iy, axis = 1)
val = np.sum(Xtheta * muhat_iy) - np.sum(max_i) - sigma * np.sum(np.log(d_i))
return -val
def grad_log_likelihood(theta):
nbK = np.asarray(theta).shape[0]
Xtheta = Phi_iy_k.dot(theta)/sigma
Xthetamat_iy = Xtheta.reshape(nbchoices, ninds).T
max_i = np.amax(Xthetamat_iy, axis = 1)
expPhi_iy = np.exp((Xthetamat_iy.T -max_i).T)
d_i = np.sum(expPhi_iy, axis = 1)
temp_mat = (Phi_iy_k.T * expPhi_iy.T.flatten()).T
list_temp = []
for i in range(nbchoices):
list_temp.append(temp_mat[i*ninds:(i+1)*ninds,])
n_i_k = np.sum(list_temp,axis = 0)
thegrad = muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten() - np.sum(n_i_k.T/d_i, axis = 1)
return -thegrad
theta0 = np.repeat(0,nbK)
sigma = 1
outcome = optimize.minimize(log_likelihood,method = 'CG',jac = grad_log_likelihood, x0 = theta0)
outcome
fun: 280.9103677250631 jac: array([-1.56394686e-06, -3.55409114e-07, -1.80343209e-06, 1.80183436e-07, 4.66669048e-07, 2.80306917e-07]) message: 'Optimization terminated successfully.' nfev: 31 nit: 16 njev: 31 status: 0 success: True x: array([-0.02769993, -0.35847009, 0.01487751, 0.01416175, 0.11784738, -0.25344193])
temp_mle = 1 / outcome['x'][nbK - 1]
theta_mle = outcome['x']*temp_mle
print(temp_mle)
print(theta_mle)
-3.9456769581766116 [ 0.10929498 1.41440717 -0.05870185 -0.05587768 -0.46498768 1. ]
lenobj = nbK+ninds
c = np.concatenate((muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten(),np.repeat(-1,ninds)))
m = grb.Model('lp')
x = m.addMVar(lenobj, name='x', lb=-grb.GRB.INFINITY)
m.setObjective(c @ x, grb.GRB.MAXIMIZE)
cstMat = spr.hstack((spr.csr_matrix(-Phi_iy_k), spr.kron(two_d(np.ones(nbchoices)),spr.identity(ninds))))
rhs = np.zeros(ninds*nbchoices)
m.addConstr(cstMat @ x >= rhs)
nbCstr = cstMat.shape[0]
const_2 = np.array([0]*(nbK - 1))
const_2 = np.append(const_2, 1)
const_2 = np.append(const_2 ,[0]*ninds)
m.addConstr(const_2 @ x == 1)
m.optimize()
if m.status == grb.GRB.Status.OPTIMAL:
print("Value of the problem (Gurobi) =", m.objval)
opt_x = m.getAttr('x')
Academic license - for non-commercial use only - expires 2022-11-14 Using license file c:\gurobi911\gurobi.lic Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64) Thread count: 6 physical cores, 12 logical processors, using up to 12 threads Optimize a model with 841 rows, 216 columns and 5881 nonzeros Model fingerprint: 0xece1c49f Coefficient statistics: Matrix range [3e-03, 5e+00] Objective range [2e-01, 5e+01] Bounds range [0e+00, 0e+00] RHS range [1e+00, 1e+00] Presolve removed 211 rows and 1 columns Presolve time: 0.01s Presolved: 630 rows, 215 columns, 2520 nonzeros Iteration Objective Primal Inf. Dual Inf. Time 0 handle free variables 0s 380 -1.3829488e+02 0.000000e+00 0.000000e+00 0s Solved in 380 iterations and 0.02 seconds Optimal objective -1.382948769e+02 Value of the problem (Gurobi) = -138.29487687593996
theta_lp = np.array(opt_x[:nbK])
print(theta_lp)
print(theta_mle)
[ 0.13949426 0.05032141 0.03395604 -0.72749334 -0.02506117 1. ] [ 0.10929498 1.41440717 -0.05870185 -0.05587768 -0.46498768 1. ]
indMax=100
tempMax=temp_mle
outcomemat = np.zeros((indMax+1,nbK-1))
def log_likelihood_fixedtemp(subsetoftheta, *temp):
val = log_likelihood(np.append(subsetoftheta, 1/temp[0]))
return val
def grad_log_likelihood_fixedtemp(subsetoftheta, *temp):
val = grad_log_likelihood(np.append(subsetoftheta, 1/temp[0]))[:-1]
return val
outcomemat[0,:] = theta_lp[:-1]
iterMax = indMax+1
for k in range(2,iterMax+1,1):
thetemp = tempMax * (k-1)/indMax
outcomeFixedTemp = optimize.minimize(log_likelihood_fixedtemp,method = 'CG',jac = grad_log_likelihood_fixedtemp, args = (thetemp,), x0 = theta0[:-1])
outcomemat[k-1,:] = outcomeFixedTemp['x']*thetemp
outcomemat
array([[ 0.13949426, 0.05032141, 0.03395604, -0.72749334, -0.02506117], [ 0.14292345, 0.42333439, -0.06371649, -0.51242318, -0.1314772 ], [ 0.14482454, 0.42218714, -0.06321331, -0.50753462, -0.12936482], [ 0.14587841, 0.42033993, -0.06257653, -0.5042307 , -0.13058508], [ 0.14624026, 0.41848991, -0.06210426, -0.50100175, -0.13399352], [ 0.14642665, 0.41756089, -0.06139318, -0.49794246, -0.13827628], [ 0.14665808, 0.41780041, -0.06035416, -0.49509681, -0.14277605], [ 0.14695082, 0.41912092, -0.05907464, -0.4924095 , -0.14720936], [ 0.1472558 , 0.42138082, -0.05767761, -0.48979281, -0.15145913], [ 0.1475228 , 0.42447024, -0.05626612, -0.48716609, -0.15548361], [ 0.14771955, 0.42831363, -0.05491026, -0.48446871, -0.15928106], [ 0.147832 , 0.43285517, -0.05365016, -0.48166091, -0.16287184], [ 0.14785899, 0.43804697, -0.05250383, -0.47872002, -0.16628727], [ 0.14780664, 0.44384285, -0.05147502, -0.4756363 , -0.16956138], [ 0.14768441, 0.4501967 , -0.05055895, -0.47240827, -0.1727272 ], [ 0.14750252, 0.45706218, -0.049747 , -0.46904023, -0.17581347], [ 0.14727078, 0.46439392, -0.04902871, -0.46553968, -0.17884437], [ 0.14699769, 0.47214821, -0.04839391, -0.46191531, -0.18183989], [ 0.14669068, 0.48028388, -0.04783308, -0.45817657, -0.18481554], [ 0.14635581, 0.4887628 , -0.04733761, -0.454333 , -0.18778329], [ 0.14599821, 0.49755006, -0.04690026, -0.45039344, -0.19075232], [ 0.14562192, 0.50661408, -0.04651459, -0.44636624, -0.19372942], [ 0.14523021, 0.51592631, -0.04617522, -0.44225927, -0.1967194 ], [ 0.14482593, 0.52546154, -0.04587746, -0.4380792 , -0.19972597], [ 0.14441128, 0.53519698, -0.04561735, -0.43383239, -0.20275147], [ 0.14398796, 0.54511243, -0.04539156, -0.4295247 , -0.2057973 ], [ 0.1435575 , 0.55519007, -0.0451971 , -0.4251611 , -0.20886437], [ 0.1431211 , 0.56541397, -0.0450314 , -0.42074623, -0.21195309], [ 0.14267966, 0.57577012, -0.0448923 , -0.41628419, -0.21506341], [ 0.14223423, 0.58624609, -0.04477745, -0.41177899, -0.21819488], [ 0.14178533, 0.59683058, -0.04468552, -0.40723351, -0.22134761], [ 0.14133354, 0.60751392, -0.04461472, -0.40265106, -0.2245206 ], [ 0.14087921, 0.61828703, -0.04456364, -0.3980343 , -0.22771335], [ 0.14042282, 0.62914224, -0.0445311 , -0.3933857 , -0.23092522], [ 0.13996489, 0.64007282, -0.04451558, -0.38870766, -0.23415508], [ 0.13950548, 0.65107205, -0.04451636, -0.38400201, -0.23740251], [ 0.13904476, 0.66213441, -0.0445324 , -0.37927088, -0.2406666 ], [ 0.13858306, 0.67325491, -0.04456273, -0.37451579, -0.24394668], [ 0.13812065, 0.68442903, -0.04460643, -0.36973808, -0.24724221], [ 0.13765741, 0.69565258, -0.04466292, -0.36493967, -0.25055202], [ 0.13719355, 0.70692191, -0.04473144, -0.36012167, -0.25387558], [ 0.13672926, 0.71823353, -0.04481122, -0.35528524, -0.25721237], [ 0.13626446, 0.7295844 , -0.04490188, -0.35043145, -0.26056178], [ 0.13579942, 0.74097186, -0.04500259, -0.34556132, -0.26392307], [ 0.13533403, 0.7523934 , -0.04511293, -0.34067596, -0.26729566], [ 0.13486843, 0.76384653, -0.04523245, -0.33577616, -0.27067909], [ 0.13440265, 0.77532938, -0.04536059, -0.3308628 , -0.27407258], [ 0.13393689, 0.78684001, -0.04549687, -0.32593626, -0.27747614], [ 0.13347077, 0.7983763 , -0.04564119, -0.32099779, -0.28088885], [ 0.1330048 , 0.809937 , -0.04579276, -0.31604763, -0.28431068], [ 0.13253869, 0.82152068, -0.04595141, -0.31108663, -0.28774075], [ 0.13207258, 0.83312551, -0.04611692, -0.30611517, -0.2911789 ], [ 0.13160643, 0.84475064, -0.04628883, -0.30113389, -0.29462481], [ 0.13114011, 0.85639463, -0.04646697, -0.29614321, -0.29807827], [ 0.13067406, 0.86805665, -0.04665084, -0.29114378, -0.30153838], [ 0.13020788, 0.87973545, -0.04684038, -0.28613587, -0.30500532], [ 0.12974168, 0.89143024, -0.04703527, -0.28112007, -0.30847847], [ 0.12927578, 0.90314039, -0.04723521, -0.2760961 , -0.3119581 ], [ 0.12880981, 0.91486466, -0.04744003, -0.27106525, -0.31544328], [ 0.12834383, 0.92660242, -0.04764964, -0.26602702, -0.31893441], [ 0.12787795, 0.93835309, -0.04786373, -0.26098232, -0.32243073], [ 0.1274121 , 0.95011595, -0.04808207, -0.25593121, -0.3259322 ], [ 0.12694638, 0.9618904 , -0.04830442, -0.25087406, -0.32943869], [ 0.12648065, 0.97367593, -0.04853098, -0.24581086, -0.3329498 ], [ 0.12601492, 0.98547194, -0.04876116, -0.2407422 , -0.33646546], [ 0.12554925, 0.99727762, -0.04899523, -0.23566864, -0.33998514], [ 0.12508389, 1.00909342, -0.04923243, -0.230589 , -0.34350961], [ 0.12461853, 1.02091806, -0.0494732 , -0.22550493, -0.34703787], [ 0.12415302, 1.03275119, -0.04971724, -0.22041607, -0.35057007], [ 0.12368792, 1.04459299, -0.04996427, -0.21532268, -0.35410574], [ 0.12322275, 1.05644253, -0.05021434, -0.2102249 , -0.35764511], [ 0.12275764, 1.06829964, -0.05046731, -0.20512296, -0.3611879 ], [ 0.1222923 , 1.08016371, -0.05072307, -0.20001745, -0.36473351], [ 0.12182733, 1.092035 , -0.05098155, -0.19490767, -0.36828271], [ 0.12136239, 1.10391293, -0.05124256, -0.18979413, -0.37183499], [ 0.12089771, 1.11579742, -0.05150596, -0.18467661, -0.37539059], [ 0.12043258, 1.12768748, -0.05177211, -0.1795563 , -0.37894824], [ 0.11996796, 1.13958399, -0.05204018, -0.17443247, -0.38250907], [ 0.11950335, 1.15148596, -0.05231066, -0.16930524, -0.38607249], [ 0.11903875, 1.16339335, -0.05258328, -0.16417497, -0.3896385 ], [ 0.1185742 , 1.17530596, -0.05285799, -0.15904171, -0.39320697], [ 0.11810972, 1.18722358, -0.05313473, -0.15390554, -0.39677782], [ 0.1176453 , 1.19914604, -0.05341341, -0.14876658, -0.40035098], [ 0.11718094, 1.21107314, -0.05369399, -0.14362492, -0.40392636], [ 0.11671664, 1.22300473, -0.05397639, -0.13848065, -0.40750389], [ 0.1162524 , 1.23494063, -0.05426056, -0.13333385, -0.41108351], [ 0.11578821, 1.2468807 , -0.05454644, -0.12818461, -0.41466514], [ 0.11532409, 1.25882477, -0.05483398, -0.123033 , -0.41824872], [ 0.11486001, 1.27077273, -0.05512312, -0.1178791 , -0.42183419], [ 0.114396 , 1.28272443, -0.05541383, -0.11272299, -0.42542148], [ 0.11393203, 1.29467973, -0.05570604, -0.10756472, -0.42901055], [ 0.11346804, 1.30663857, -0.05599965, -0.10240454, -0.43260135], [ 0.11300418, 1.31860074, -0.05629474, -0.0972422 , -0.4361938 ], [ 0.11254036, 1.33056618, -0.05659122, -0.09207789, -0.43978786], [ 0.11207659, 1.34253477, -0.05688903, -0.08691169, -0.44338348], [ 0.11161287, 1.35450642, -0.05718813, -0.08174365, -0.44698063], [ 0.11114919, 1.36648102, -0.0574885 , -0.07657383, -0.45057925], [ 0.11068556, 1.37845848, -0.05779009, -0.07140228, -0.4541793 ], [ 0.11022198, 1.39043872, -0.05809287, -0.06622904, -0.45778073], [ 0.10975844, 1.40242164, -0.05839681, -0.06105416, -0.46138352], [ 0.10929494, 1.41440717, -0.05870187, -0.0558777 , -0.46498762]])
The zero-temperature estimator is:
print(outcomemat[1,:])
[ 0.14292345 0.42333439 -0.06371649 -0.51242318 -0.1314772 ]
The mle estimator is:
print(outcomemat[indMax,])
[ 0.10929494 1.41440717 -0.05870187 -0.0558777 -0.46498762]
nbB = 2
thetemp = 1
epsilon_biy = special.digamma(1) -np.log(-np.log(np.random.uniform(0,1,ninds*nbchoices*nbB)))
lenobj = ninds*nbB+nbK
newc = np.concatenate((muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten(),np.repeat(-1/nbB,ninds*nbB)))
newm = grb.Model('new_lp')
x = newm.addMVar(lenobj, name='x', lb=-grb.GRB.INFINITY)
newm.setObjective(newc @ x, grb.GRB.MAXIMIZE)
mat1 = spr.kron(-Phi_iy_k, two_d(np.repeat(1,nbB)))
mat2 = spr.kron(two_d(np.repeat(1,nbchoices)),spr.identity(ninds*nbB))
newcstMat = spr.hstack((mat1, mat2))
rhs = epsilon_biy
newm.addConstr(newcstMat @ x >= rhs)
newm.optimize()
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64) Thread count: 6 physical cores, 12 logical processors, using up to 12 threads Optimize a model with 1680 rows, 426 columns and 11760 nonzeros Model fingerprint: 0x93e351cb Coefficient statistics: Matrix range [3e-03, 5e+00] Objective range [2e-01, 5e+01] Bounds range [0e+00, 0e+00] RHS range [5e-04, 7e+00] Presolve time: 0.01s Presolved: 1680 rows, 426 columns, 11760 nonzeros Iteration Objective Primal Inf. Dual Inf. Time 0 handle free variables 0s 1123 -2.9966123e+02 0.000000e+00 0.000000e+00 0s Solved in 1123 iterations and 0.06 seconds Optimal objective -2.996612347e+02
if m.status == grb.GRB.Status.OPTIMAL:
print("Value of the problem (Gurobi) =", newm.objval)
opt_x = np.array(newm.getAttr('x'))
newtheta_lp = opt_x[:nbK] / opt_x[nbK-1]
Value of the problem (Gurobi) = -299.661234703907
print(theta_mle)
print(newtheta_lp)
[ 0.10929498 1.41440717 -0.05870185 -0.05587768 -0.46498768 1. ] [-0.01345289 1.2917597 -0.22401021 -0.00305825 -0.50896292 1. ]
Finally probit
nbB = 2
thetemp = 1
epsilon_biy = np.random.normal(nbB*ninds*nbchoices)
lenobj = ninds*nbB+nbK
newc = np.concatenate((muhat_iy.reshape(1,nbchoices*ninds).dot(Phi_iy_k).flatten(),np.repeat(-1/nbB,ninds*nbB)))
newm = grb.Model('new_lp')
x = newm.addMVar(lenobj, name='x', lb=-grb.GRB.INFINITY)
newm.setObjective(newc @ x, grb.GRB.MAXIMIZE)
mat1 = spr.kron(-Phi_iy_k, two_d(np.repeat(1,nbB)))
mat2 = spr.kron(two_d(np.repeat(1,nbchoices)),spr.identity(ninds*nbB))
newcstMat = spr.hstack((mat1, mat2))
rhs = epsilon_biy
newm.addConstr(newcstMat @ x >= rhs)
newm.optimize()
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64) Thread count: 6 physical cores, 12 logical processors, using up to 12 threads Optimize a model with 1680 rows, 426 columns and 11760 nonzeros Model fingerprint: 0x5e8fff22 Coefficient statistics: Matrix range [3e-03, 5e+00] Objective range [2e-01, 5e+01] Bounds range [0e+00, 0e+00] RHS range [2e+03, 2e+03] Presolve time: 0.01s Presolved: 1680 rows, 426 columns, 11760 nonzeros Iteration Objective Primal Inf. Dual Inf. Time 0 handle free variables 0s 750 -3.5266815e+05 0.000000e+00 0.000000e+00 0s Solved in 750 iterations and 0.04 seconds Optimal objective -3.526681518e+05