# %% [markdown] # # Discrete choice experiment: preferences for a digital savings product # # **Worked example on synthetic data, by Ahmed Beddi.** The respondents and their choices are simulated # from tastes chosen for illustration; no survey was fielded and no real product, provider or market is # represented. Money amounts are in generic local currency units (LCU). # # Contents: (1) attributes and levels; (2) an efficient design found by coordinate exchange on the # D-error of the multinomial logit, with blocking; (3) simulated respondents with heterogeneous tastes; # (4) conditional logit by maximum likelihood (analytic gradient and Hessian); (5) a mixed logit sketch # by simulated maximum likelihood; (6) willingness to pay with delta-method and Krinsky–Robb intervals; # (7) market-share simulation, including the numbers used by the website demo. # %% import itertools, math import numpy as np import pandas as pd from scipy.optimize import minimize from scipy.stats import norm, qmc pd.set_option('display.width', 120) rng = np.random.default_rng(20261103) # %% [markdown] # ## 1. Attributes and levels # # | Attribute | Levels | Coding | # |---|---|---| # | Monthly fee | 0, 100, 200, 400 LCU | continuous, per 100 LCU | # | Withdrawal rule | flexible; early withdrawal with a penalty ("soft lock"); locked until the goal date ("strict lock") | dummies, base = flexible | # | Walk to the nearest cash-out agent | 5, 20, 45 minutes | continuous, per 10 minutes | # | Reminder | none; SMS text; voice message | dummies, base = none | # | Bonus at goal date | 0, 3, 6 % of balance | continuous, per percentage point | # # Each choice task shows two unlabelled products (A, B) and an opt-out: *keep my current way of saving*. # The opt-out carries an alternative-specific constant (ASC). Without it, shares cannot be predicted, # because people are forced into a product they may not want. # %% LEVELS = { 'fee': [0, 100, 200, 400], 'withdraw': [0, 1, 2], # 0 flexible, 1 soft lock, 2 strict lock 'distance': [5, 20, 45], 'reminder': [0, 1, 2], # 0 none, 1 SMS, 2 voice 'bonus': [0, 3, 6], } ATTRS = list(LEVELS) NAMES = ['fee_100', 'soft_lock', 'strict_lock', 'dist_10min', 'sms', 'voice', 'bonus_pct', 'asc_optout'] K = len(NAMES) def code_alt(a): """Level vector (fee, withdraw, distance, reminder, bonus) -> coded row (without the ASC).""" fee, w, d, r, b = a return np.array([fee / 100, w == 1, w == 2, d / 10, r == 1, r == 2, b, 0.0], float) OPTOUT = np.eye(K)[-1] def design_to_X(D): """D: (S, 2, 5) array of levels -> X: (S, 3, K) coded, third alternative = opt-out.""" X = np.zeros((D.shape[0], 3, K)) X[:, :2, 0] = D[:, :, 0] / 100 X[:, :2, 1] = D[:, :, 1] == 1 X[:, :2, 2] = D[:, :, 1] == 2 X[:, :2, 3] = D[:, :, 2] / 10 X[:, :2, 4] = D[:, :, 3] == 1 X[:, :2, 5] = D[:, :, 3] == 2 X[:, :2, 6] = D[:, :, 4] X[:, 2] = OPTOUT return X # %% [markdown] # ## 2. Efficient design: D-error and coordinate exchange # # For the multinomial logit with prior parameters $\tilde\beta$, the Fisher information of a design with # choice sets $s$ and alternatives $j$ is # # $$I(\tilde\beta) = \sum_s \sum_j P_{sj}\,(x_{sj} - \bar x_s)(x_{sj} - \bar x_s)^\top,\qquad \bar x_s = \sum_j P_{sj}\,x_{sj},$$ # # and the D-error is $\det\big(I(\tilde\beta)^{-1}\big)^{1/K}$: smaller is better (smaller joint confidence # region for the parameters). Priors are deliberately modest and only sign-informed, as one would set them # before a pilot: a design optimised for confident but wrong priors can be worse than one with small priors. # # Two searches, both rejecting any change that makes one product dominate the other (at least as good on # fee, distance and bonus, identical on the rest): # # * **Coordinate exchange** (unconstrained): for each choice set, alternative and attribute in turn, try # every level and keep the one that lowers the D-error. For attributes coded linearly this piles the # design onto the extreme levels, which is efficient only if effects are truly linear. # * **Swapping** (Huber and Zwerina, 1996): start from a random design in which every level appears equally # often, and only *swap* levels of one attribute between two alternatives. Level balance is preserved by # construction, so middle levels stay in the design and non-linearity can be checked. This is the design # used below. # # Each search restarts from 20 random designs and keeps the best. # %% PRIOR = np.array([-0.25, 0.0, -0.1, -0.1, 0.15, 0.2, 0.1, 0.5]) S, BLOCKS = 18, 3 def info(X, beta): U = X @ beta P = np.exp(U - U.max(1, keepdims=True)); P /= P.sum(1, keepdims=True) xbar = (P[:, :, None] * X).sum(1) Z = X - xbar[:, None, :] return np.einsum('sj,sjk,sjl->kl', P, Z, Z) def d_error(D, beta=PRIOR): I = info(design_to_X(D), beta) sign, logdet = np.linalg.slogdet(I) return np.inf if sign <= 0 else math.exp(-logdet / K) def dominated(a, b): """True if one alternative is at least as good on fee, distance and bonus, strictly better on one, and identical on the unsigned attributes.""" if a[1] != b[1] or a[3] != b[3]: return False ge = (a[0] <= b[0]) and (a[2] <= b[2]) and (a[4] >= b[4]) le = (a[0] >= b[0]) and (a[2] >= b[2]) and (a[4] <= b[4]) return (ge or le) and tuple(a) != tuple(b) or tuple(a) == tuple(b) def random_design(g): D = np.zeros((S, 2, 5), int) for j in range(2): for k, a in enumerate(ATTRS): lv = np.resize(LEVELS[a], S); D[:, j, k] = g.permutation(lv) return D def coordinate_exchange(g, max_iter=20): D = random_design(g) while any(dominated(D[s, 0], D[s, 1]) for s in range(S)): D = random_design(g) best = d_error(D) for _ in range(max_iter): improved = False for s in range(S): for j in range(2): for k, a in enumerate(ATTRS): cur = D[s, j, k] for lv in LEVELS[a]: if lv == cur: continue D[s, j, k] = lv if dominated(D[s, 0], D[s, 1]): D[s, j, k] = cur; continue e = d_error(D) if e < best - 1e-12: best, cur, improved = e, lv, True else: D[s, j, k] = cur D[s, j, k] = cur if not improved: break return D, best def swap_search(g, passes=8): D = random_design(g) while any(dominated(D[s, 0], D[s, 1]) for s in range(S)): D = random_design(g) best = d_error(D) cells = [(s, j) for s in range(S) for j in range(2)] for _ in range(passes): improved = False for k in range(len(ATTRS)): for a in range(len(cells)): for b in range(a + 1, len(cells)): (s1, j1), (s2, j2) = cells[a], cells[b] if D[s1, j1, k] == D[s2, j2, k]: continue D[s1, j1, k], D[s2, j2, k] = D[s2, j2, k], D[s1, j1, k] if dominated(D[s1, 0], D[s1, 1]) or dominated(D[s2, 0], D[s2, 1]): e = np.inf else: e = d_error(D) if e < best - 1e-12: best, improved = e, True else: D[s1, j1, k], D[s2, j2, k] = D[s2, j2, k], D[s1, j1, k] if not improved: break return D, best g = np.random.default_rng(7) D_ce, derr_ce = min((coordinate_exchange(g) for _ in range(20)), key=lambda t: t[1]) D, derr = min((swap_search(g) for _ in range(20)), key=lambda t: t[1]) rand_err = np.array([d_error(random_design(g)) for _ in range(500)]) rand_err = rand_err[np.isfinite(rand_err)] print(f'random level-balanced designs: median D-error {np.median(rand_err):.4f} (best of 500: {rand_err.min():.4f})') print(f'swapping, level-balanced (chosen): D-error {derr:.4f} D-efficiency vs median random: {np.median(rand_err) / derr:.2f}') print(f'coordinate exchange, unconstrained: D-error {derr_ce:.4f} (extreme levels only, see below)') # %% [markdown] # Level frequencies: the unconstrained design drops the middle levels of the linear attributes; the # swapped design keeps every level. # %% for lab, DD in [('coordinate exchange', D_ce), ('swapping (chosen)', D)]: print(lab) for k, a in enumerate(ATTRS): vals, cnt = np.unique(DD[:, :, k], return_counts=True) print(f' {a:<9}', dict(zip(vals.tolist(), cnt.tolist()))) # %% [markdown] # ### Blocking # # Eighteen tasks are too many for respondents with limited schooling answering by voice or with picture # cards; each respondent sees six. The 18 sets are split into 3 blocks of 6 so that, within each block, # every attribute level appears as close as possible to its overall share (random search over # 20,000 partitions). # %% def block_score(assign): sc = 0.0 for k, a in enumerate(ATTRS): for lv in LEVELS[a]: tot = (D[:, :, k] == lv).sum() / BLOCKS for b in range(BLOCKS): sc += ((D[assign == b, :, k] == lv).sum() - tot) ** 2 return sc gb = np.random.default_rng(11) base = np.repeat(np.arange(BLOCKS), S // BLOCKS) best_b, best_sc = None, np.inf for _ in range(20000): a = gb.permutation(base); sc = block_score(a) if sc < best_sc: best_b, best_sc = a, sc block = best_b print('block score (sum of squared deviations):', round(best_sc, 2)) rows = [] lab_w = ['flexible', 'soft lock', 'strict lock']; lab_r = ['none', 'SMS', 'voice'] order = np.lexsort((np.arange(S), block)) for new_id, s in enumerate(order, start=1): for j, alt in enumerate('AB'): f, w, d, r, b = D[s, j] rows.append(dict(set=new_id, block=int(block[s]) + 1, alt=alt, fee=int(f), withdrawal=lab_w[w], walk_min=int(d), reminder=lab_r[r], bonus_pct=int(b))) design = pd.DataFrame(rows) design.to_csv('dce_design.csv', index=False) D, block = D[order], block[order] print(design.head(12).to_string(index=False)) # %% [markdown] # ## 3. Simulated respondents # # 600 respondents, 200 per block, half with low literacy. True tastes (mixed logit, chosen for # illustration): # # * fee −0.35 per 100 LCU; distance −0.12 per 10 minutes; bonus +0.15 per percentage point (fixed); # * soft lock $\sim N(0.15, 0.5^2)$; strict lock $\sim N(-0.30, 1.0^2)$: commitment is attractive to some # people and repellent to others, as commitment-savings evidence suggests # (Ashraf, Karlan and Yin, 2006); # * SMS +0.35 if literate, +0.05 if low literacy; voice +0.20 if literate, +0.60 if low literacy; # * opt-out constant +0.8. # # Errors are i.i.d. type-I extreme value, so choices follow a logit kernel. The same function generates # the replications in section 6. # %% N_RESP, TASKS = 600, 6 TRUE = dict(fee=-0.35, dist=-0.12, bonus=0.15, soft_m=0.15, soft_s=0.5, strict_m=-0.30, strict_s=1.0, sms_lit=0.35, sms_low=0.05, voice_lit=0.20, voice_low=0.60, asc=0.8) Xd = design_to_X(D) def simulate_choices(g): low = g.random(N_RESP) < 0.5 blk = np.repeat([0, 1, 2], N_RESP // 3) soft = g.normal(TRUE['soft_m'], TRUE['soft_s'], N_RESP); strict = g.normal(TRUE['strict_m'], TRUE['strict_s'], N_RESP) B = np.column_stack([np.full(N_RESP, TRUE['fee']), soft, strict, np.full(N_RESP, TRUE['dist']), np.where(low, TRUE['sms_low'], TRUE['sms_lit']), np.where(low, TRUE['voice_low'], TRUE['voice_lit']), np.full(N_RESP, TRUE['bonus']), np.full(N_RESP, TRUE['asc'])]) sets = np.array([np.flatnonzero(block == b) for b in range(3)])[blk] # (N, 6) U = np.einsum('nsjk,nk->nsj', Xd[sets], B) + g.gumbel(size=(N_RESP, TASKS, 3)) df = pd.DataFrame({'id': np.repeat(np.arange(N_RESP), TASKS), 'lowlit': np.repeat(low.astype(int), TASKS), 'set': sets.ravel(), 'choice': U.argmax(2).ravel()}) return df, strict ch, strict_draws = simulate_choices(rng) print(f'realised mean of strict-lock tastes in this sample: {np.mean(strict_draws):.3f} (population mean -0.30)') print('observations:', len(ch), ' choice shares (A, B, opt-out):', ch.choice.value_counts(normalize=True).sort_index().round(3).tolist()) # %% [markdown] # ## 4. Conditional logit by maximum likelihood # # $$P_{nsj} = \frac{\exp(x_{sj}^\top\beta)}{\sum_{l} \exp(x_{sl}^\top\beta)},\qquad # \ell(\beta) = \sum_{n,s} \log P_{ns,y_{ns}},\qquad # \nabla\ell = \sum_{n,s}\Big(x_{s y_{ns}} - \sum_j P_{nsj} x_{sj}\Big).$$ # # The Hessian is $-\sum_{n,s}\sum_j P_{nsj}(x_{sj}-\bar x_{ns})(x_{sj}-\bar x_{ns})^\top$; standard errors # come from its inverse. Because each respondent answers six tasks, and her tastes persist across them, # the reported standard errors are sandwich estimates clustered by respondent; the naive inverse-Hessian # errors are shown for comparison and are too small. Model (a) has main effects only; model (b) adds SMS × low-literacy and # voice × low-literacy interactions. The demo uses model (b). # %% def build(ch, interactions): X = Xd[ch.set.values].copy() if interactions: low = ch.lowlit.values[:, None].astype(float) X = np.concatenate([X, (X[:, :, 4] * low)[:, :, None], (X[:, :, 5] * low)[:, :, None]], axis=2) return X, ch.choice.values def cl_fit(X, y, groups=None): def nll(b): U = X @ b; U -= U.max(1, keepdims=True) lse = np.log(np.exp(U).sum(1)) return -(U[np.arange(len(y)), y] - lse).sum() def grad(b): U = X @ b; P = np.exp(U - U.max(1, keepdims=True)); P /= P.sum(1, keepdims=True) return -(X[np.arange(len(y)), y] - (P[:, :, None] * X).sum(1)).sum(0) r = minimize(nll, np.zeros(X.shape[2]), jac=grad, method='BFGS', options={'gtol': 1e-8}) b = r.x U = X @ b; P = np.exp(U - U.max(1, keepdims=True)); P /= P.sum(1, keepdims=True) Z = X - (P[:, :, None] * X).sum(1)[:, None, :] H = np.einsum('nj,njk,njl->kl', P, Z, Z) V = np.linalg.inv(H) if groups is not None: # sandwich clustered by respondent: tasks of one person are not independent sc = X[np.arange(len(y)), y] - (P[:, :, None] * X).sum(1) G = np.zeros((groups.max() + 1, X.shape[2])); np.add.at(G, groups, sc) V = V @ (G.T @ G) @ V return b, V, -r.fun Xa, y = build(ch, False) ba, Va, lla = cl_fit(Xa, y, ch.id.values) Xb, _ = build(ch, True) bb, Vb0, llb = cl_fit(Xb, y) _, Vb, _ = cl_fit(Xb, y, ch.id.values) names_b = NAMES + ['sms_x_lowlit', 'voice_x_lowlit'] tab = pd.DataFrame({'(a) coef': np.r_[ba, [np.nan, np.nan]], '(a) se': np.r_[np.sqrt(np.diag(Va)), [np.nan, np.nan]], '(b) coef': bb, '(b) se': np.sqrt(np.diag(Vb)), '(b) naive se': np.sqrt(np.diag(Vb0))}, index=names_b) print(tab.round(3)) print(f'log-likelihood (a) {lla:.1f} (b) {llb:.1f} LR test of interactions: chi2(2) = {2 * (llb - lla):.1f}') # %% [markdown] # Coefficients for the website demo (model b), printed in the form hardcoded in `DceExplorer.tsx`: # %% print('const BETA = {' + ', '.join(f'{n}: {v:.4f}' for n, v in zip( ['fee', 'soft', 'strict', 'dist', 'sms', 'voice', 'bonus', 'asc', 'smsLow', 'voiceLow'], bb)) + ' };') # %% [markdown] # ## 5. Mixed logit sketch (simulated maximum likelihood) # # The conditional logit gives the *average* taste for a strict lock. The simulated data were generated # with a wide spread around that average. A mixed logit lets the soft- and strict-lock coefficients be # normal across respondents, $\beta_n^{\text{lock}} \sim N(\mu, \sigma^2)$, and integrates over that # distribution respondent by respondent (panel structure: the same person keeps the same taste across her # six tasks): # # $$L_n(\theta) = \int \prod_{s} P_{ns y_{ns}}(\beta)\, \phi(\beta\mid\mu,\sigma)\,d\beta # \approx \frac{1}{R}\sum_{r=1}^{R}\prod_s P_{ns y_{ns}}(\beta^{(r)}_n),$$ # # with $R = 200$ scrambled Halton draws per respondent (Train, 2009, ch. 6 and 9). # %% ids = ch.id.values; uid, inv = np.unique(ids, return_inverse=True) nR, R = len(uid), 200 H = qmc.Halton(d=2, scramble=True, seed=3).random(nR * R).reshape(nR, R, 2) E = norm.ppf(H) # standard normal draws, (respondent, draw, 2) def mxl_nll(theta): b = theta[:10].copy(); s_soft, s_strict = np.abs(theta[10]), np.abs(theta[11]) Ufix = Xb @ b # (obs, 3) eps_soft = E[inv, :, 0] * s_soft # (obs, R) eps_strict = E[inv, :, 1] * s_strict U = Ufix[:, None, :] + eps_soft[:, :, None] * Xb[:, None, :, 1] + eps_strict[:, :, None] * Xb[:, None, :, 2] U -= U.max(2, keepdims=True) logP = U[np.arange(len(y)), :, y] - np.log(np.exp(U).sum(2)) # (obs, R) lp = np.zeros((nR, R)); np.add.at(lp, inv, logP) # product over a person's tasks m = lp.max(1, keepdims=True) return -(m[:, 0] + np.log(np.exp(lp - m).mean(1))).sum() th0 = np.r_[bb, 0.3, 0.3] rm = minimize(mxl_nll, th0, method='L-BFGS-B', options={'maxiter': 500}) th = rm.x.copy(); th[10:] = np.abs(th[10:]) # numerical Hessian for standard errors def num_hess(f, x, h=1e-4): k = len(x); Hm = np.zeros((k, k)) for i in range(k): for j in range(i, k): e_i, e_j = np.eye(k)[i] * h, np.eye(k)[j] * h Hm[i, j] = Hm[j, i] = (f(x + e_i + e_j) - f(x + e_i - e_j) - f(x - e_i + e_j) + f(x - e_i - e_j)) / (4 * h * h) return Hm Vm = np.linalg.inv(num_hess(mxl_nll, th)) mxl = pd.DataFrame({'coef': th, 'se': np.sqrt(np.diag(Vm))}, index=names_b + ['sd_soft_lock', 'sd_strict_lock']) print(mxl.round(3)) print(f'simulated log-likelihood {-rm.fun:.1f} vs conditional logit (b) {llb:.1f}') print('true values used in simulation: soft lock N(0.15, 0.5), strict lock N(-0.30, 1.0)') # %% [markdown] # The mixed logit recovers a large standard deviation for the strict lock: a sizeable minority values # a locked product positively even though the average taste is negative, and the conditional logit # average hides this. Mixed-logit means are also on a different scale from conditional-logit coefficients # (the unobserved taste variation is no longer in the error term), so their absolute sizes are not # directly comparable; ratios such as WTP are. # # ## 6. Willingness to pay # # WTP for attribute $k$ in LCU per month is $\text{WTP}_k = -100\,\beta_k/\beta_{\text{fee}}$ (the fee is # coded per 100 LCU). Two intervals: # # * **delta method:** $\operatorname{Var}(\widehat{\text{WTP}}_k) \approx \nabla g^\top V \nabla g$ with # $\partial g/\partial\beta_k = -100/\beta_{\text{fee}}$ and # $\partial g/\partial\beta_{\text{fee}} = 100\,\beta_k/\beta_{\text{fee}}^2$; # * **Krinsky–Robb:** draw 10,000 parameter vectors from $N(\hat\beta, \hat V)$, compute the ratio for each, # take the 2.5th and 97.5th percentiles (Krinsky and Robb, 1986). # %% def wtp_table(b, V, names, fee_idx=0, skip=('fee_100', 'asc_optout', 'sms_x_lowlit', 'voice_x_lowlit')): out = []; kr = np.random.default_rng(5).multivariate_normal(b, V, 10000) for k, n in enumerate(names): if n in skip: continue w = -100 * b[k] / b[fee_idx] gvec = np.zeros(len(b)); gvec[k] = -100 / b[fee_idx]; gvec[fee_idx] = 100 * b[k] / b[fee_idx] ** 2 se = math.sqrt(gvec @ V @ gvec) draws = -100 * kr[:, k] / kr[:, fee_idx] lo, hi = np.percentile(draws, [2.5, 97.5]) out.append(dict(attribute=n, WTP=w, delta_lo=w - 1.96 * se, delta_hi=w + 1.96 * se, KR_lo=lo, KR_hi=hi)) return pd.DataFrame(out).set_index('attribute') wt = wtp_table(bb, Vb, names_b).rename(index={'sms': 'sms (literate)', 'voice': 'voice (literate)'}) # WTP for voice and SMS among low-literacy respondents: beta_voice + beta_voice_x_lowlit L = np.zeros((2, len(bb))); L[0, 4] = L[0, 8] = 1; L[1, 5] = L[1, 9] = 1 b_low = np.r_[bb[0], L @ bb]; V_low = np.zeros((3, 3)) Tm = np.vstack([np.eye(len(bb))[0], L]); V_low = Tm @ Vb @ Tm.T wl = wtp_table(b_low, V_low, ['fee_100', 'sms (low literacy)', 'voice (low literacy)']) print(pd.concat([wt, wl]).round(0)) true_w = {'voice (literate)': -100 * 0.20 / -0.35, 'voice (low literacy)': -100 * 0.60 / -0.35, 'bonus per pct': -100 * 0.15 / -0.35} print('true WTP in the simulation:', {k: round(v) for k, v in true_w.items()}) # %% [markdown] # Reading the table. For low-literacy respondents a voice reminder is worth several percentage points of # bonus, and much more than an SMS; for literate respondents the two channels are close. The average WTP # for the strict lock is near zero and imprecise, but section 5 shows a wide spread of tastes around that # average: for a commitment feature the average is the wrong summary. The opt-out constant is not a WTP for # anything; it measures the pull of the status quo, and is left out of the table. # %% [markdown] # ### Is this one sample typical? # # One synthetic sample of 600 can land far from the truth, just as one real survey can. To see the # sampling behaviour of the estimator rather than of one draw, we repeat the whole simulation 100 times # (same design and same simulation function, new respondents and choices) and record the conditional-logit WTP and whether its 95 % # delta-method interval (respondent-clustered) covers the true value. Fee, distance and bonus are fixed # tastes in the simulation, and the voice coefficient is fixed within each literacy group, so these # ratios have a clean true value. # %% truth = {'dist_10min': -100 * TRUE['dist'] / TRUE['fee'], 'bonus_pct': -100 * TRUE['bonus'] / TRUE['fee'], 'voice (low literacy)': -100 * TRUE['voice_low'] / TRUE['fee']} mc = [] gmc = np.random.default_rng(31) for rep in range(100): c, _ = simulate_choices(gmc); Xr, yr = build(c, True) br, Vr, _ = cl_fit(Xr, yr, c.id.values) w1 = wtp_table(br, Vr, names_b) Tm_r = np.vstack([np.eye(len(br))[0], L]) w2 = wtp_table(np.r_[br[0], L @ br], Tm_r @ Vr @ Tm_r.T, ['fee_100', 'sms (low literacy)', 'voice (low literacy)']) w = pd.concat([w1, w2]) for k, tv in truth.items(): mc.append(dict(attr=k, est=w.loc[k, 'WTP'], cover=w.loc[k, 'delta_lo'] <= tv <= w.loc[k, 'delta_hi'])) mc = pd.DataFrame(mc).groupby('attr').agg(mean_WTP=('est', 'mean'), sd_WTP=('est', 'std'), coverage=('cover', 'mean')) mc['true_WTP'] = pd.Series(truth) print(mc.round(2)) # %% [markdown] # ## 7. Market-share simulation # # Shares for two configured products and the opt-out, from model (b), for each literacy segment and for # the sample mix (50 % low literacy). The defaults below are the defaults of the website demo. # %% def coded(prod, low): fee, w, d, r, b = prod x = code_alt(prod)[:7] return np.r_[x, 0.0, x[4] * low, x[5] * low] def shares(pA, pB, low, beta=bb): U = np.array([coded(pA, low) @ beta, coded(pB, low) @ beta, beta[7]]) e = np.exp(U - U.max()); return e / e.sum() A = (100, 0, 20, 1, 0) # 100 LCU, flexible, 20 min, SMS, 0 % bonus B = (200, 2, 5, 2, 6) # 200 LCU, strict lock, 5 min, voice, 6 % bonus for lab, low in [('literate', 0), ('low literacy', 1)]: print(f'{lab:<13} A / B / opt-out:', shares(A, B, low).round(3)) mix = 0.5 * shares(A, B, 0) + 0.5 * shares(A, B, 1) print('sample mix A / B / opt-out:', mix.round(3)) # %% [markdown] # **The same question with the mixed logit.** Shares are averaged over the estimated distribution of # lock tastes. Because people who dislike the strict lock now substitute toward A or the opt-out in # proportion to their own tastes (not in proportion to average shares, the IIA property of the # conditional logit), predicted shares can differ noticeably when products differ on a heterogeneous # attribute. # %% def shares_mxl(pA, pB, low, th=th, R=5000): gd = np.random.default_rng(9); b = th[:10] soft = b[1] + th[10] * gd.standard_normal(R); strict = b[2] + th[11] * gd.standard_normal(R) out = np.zeros(3) for i in range(R): bi = b.copy(); bi[1], bi[2] = soft[i], strict[i] out += shares(pA, pB, low, bi) return out / R mix_m = 0.5 * shares_mxl(A, B, 0) + 0.5 * shares_mxl(A, B, 1) print('mixed logit, sample mix A / B / opt-out:', mix_m.round(3))