<jupyter_start><jupyter_text>Simulated Probit<jupyter_code>import numpy as np
import gurobipy as grb
seed = 777
nbDraws = 10
U_y = np.array([0.4, 0.5, 0.2, 0.3, 0.1, 0])
nbY = len(U_y)
rho = 0.5
Covar = rho * np.ones((nbY, nbY)) + (1 - rho) * np.eye(nbY)
E = np.linalg.eigh(Covar)
V = E[0]
Q = E[1]
SqrtCovar = Q.dot(np.diag(np.sqrt(V))).dot(Q.T)
epsilon_iy = np.random.normal(0,1,nbDraws*nbY).reshape(nbDraws,nbY).dot(SqrtCovar)
u_iy = epsilon_iy + U_y
ui = np.max(u_iy, axis=1)
s_y = np.sum((u_iy.T - ui).T == 0, axis=0) / nbDraws
opt_assign = [(i, j) for i in range(nbY) for j in range(nbDraws)]

m = grb.Model()
vars = m.addVars(opt_assign, obj=epsilon_iy.T.ravel(), name='vars')
m.ModelSense = -1
m.addConstrs(vars.sum('*', i) == 1/nbDraws for i in range(nbDraws))
m.addConstrs(vars.sum(i, '*') == s_y[i] for i in range(nbY))
m.optimize()
if m.Status == grb.GRB.Status.OPTIMAL:
    pi = m.getAttr('pi')
    Uhat_y = -np.subtract(pi[nbDraws:nbY+nbDraws], pi[nbY + nbDraws - 1])
    print('U_y (true and recovered)')
    print(U_y)
    print(Uhat_y)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19685767097966728, 'fraction_numerical': 0.026802218114602587, 'mean_word_length': 4.3088235294117645, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '615492', 'license': 'AFL-1.1', 'n_tokens_mistral': 556, 'n_tokens_neox': 534, 'n_words': 102, 'path': 'veteran_students_mec_optim/Pauline/B12_demand.ipynb', 'repo_name': 'math-econ-code/optim_2019-06', 'size': '1082'}