Linear and logistic regression - ObjectVision/GeoDMS GitHub Wiki

Regression models can be estimated inside a GeoDMS configuration with the matrix functions: matr_var for the cross-product matrix $X'X$, matr_inv for its inverse and matr_mul for the remaining products. This page shows how, with three worked examples that run as given: a linear regression with several regressors (ordinary least squares), a binomial logistic regression and a multinomial logit model. The two logit models are estimated by Newton-Raphson steps in an iterate loop. The last section says how other discrete choice concepts, the nested logit in particular, fit the same building blocks. The theory behind the logit models is on Logit regression.

The examples were run with GeoDMS 20.19. The numbers in the tables are what the configurations produce; they agree with a numpy reference implementation to all printed digits.

A matrix is an attribute of a two-dimensional domain: a unit with an ipoint (or spoint) value type whose range has as many rows and columns as the matrix. The cells are stored row by row and the first coordinate of a point is the row, so the domain of an $N \times M$ matrix is

unit<ipoint> XM := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrRows, nrCols, ipoint));

with nrRows and nrCols int32 parameters. The literal 0i is an int32 zero; a bare 0 is a uint32 and would make a upoint. See point_yx and XY order for the coordinate order.

The three functions:

They take float32 or float64 attributes. Use float64 for regression: the cross products of a few hundred observations already exceed the seven digits of a float32.

There is no transpose function and no function that turns a set of attributes into a matrix, and neither is needed, because a two-dimensional domain and a one-dimensional domain unit with the same number of elements have the same memory layout. union_data is a concatenating relabel over any domain with the right number of elements, and a lookup through a computed point relation reaches any cell. The idioms the examples use:

• A set of attributes is the transposed data matrix. The union of the $M$ regressors of an $N$-domain, each $N$ elements long, fills an $M \times N$ matrix row by row: the first argument becomes row 0, the second row 1, and so on. So union_data(XtM, const(1.0, D), D/x1, D/x2) is $X'$, with the constant as its first row.
• Transposition is a lookup with the coordinates swapped: Xt[point_yx(pointcol(id(XM)), pointrow(id(XM)), XtM)] reads cell (column, row) of $X'$ for cell (row, column) of $X$. See pointrow and pointcol.
• A vector and a one-column matrix are relabels of each other: union_data(YM, D/y) makes the $N \times 1$ matrix, union_data(D, Yhat) brings a result back to the domain of the observations, and union_data(Coef, Bcol) turns an $M \times 1$ coefficient matrix into an attribute of the coefficient domain.
• The diagonal of an $M \times M$ matrix, for the standard errors: A[point_yx(int32(id(Coef)), int32(id(Coef)), MM)].
• A relation from the cells to the rows, to scale every row of $X$ by a weight of its observation: attribute<D> obs_rel := value(pointrow(id(.)), D) as subitem of the matrix domain, then X * sqrt(w[XM/obs_rel]).

Two remarks on names. Tree item names are case-insensitive: B and b, N and n, Bm and BM denote the same item, and the second definition is reported as already defined. A name that differs only in case from a function name or a literal suffix, Var next to the var function or B next to the uint8 suffix b, is reported as a case mix-up. The examples therefore use Coef for the coefficient domain, Bcol for the coefficient column and nrObs, nrCoef for the counts. The # prefix in #Coef is the nrofrows operator.

The linear model $y = X\beta + \varepsilon$, with $N$ observations and $M$ regressors of which the first is the constant, has the least squares estimator

$$\hat\beta = (X'X)^{-1} X'y$$

with fitted values $\hat y = X\hat\beta$, residuals $e = y - \hat y$, residual variance $s^2 = e'e / (N - M)$, covariance matrix $s^2 (X'X)^{-1}$ of the coefficients, and $R^2 = 1 - e'e / \sum_i (y_i - \bar y)^2$.

The example regresses the price of ten dwellings on their floor area and their distance to the city centre.

container regression
{
	unit<uint32> Dwelling : nrofrows = 10
	{
		attribute<float64> area     : [ 60, 75, 90, 110, 120, 85, 70, 130, 95, 105 ];      // m2
		attribute<float64> distance : [ 2.0, 5.5, 1.0, 8.0, 3.0, 6.5, 4.0, 2.5, 7.0, 1.5 ]; // km to the centre
		attribute<float64> price    : [ 220, 225, 318, 310, 392, 248, 233, 415, 284, 349 ]; // k euro
	}

	container ols
	{
		unit<uint32> Coef : nrofrows = 3 { attribute<string> name : ['const', 'area', 'distance']; }

		parameter<int32> nrObs  := int32(#Dwelling);
		parameter<int32> nrCoef := int32(#Coef);

		unit<ipoint> XM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrObs,  nrCoef, ipoint)); // N x M
		unit<ipoint> XtM := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef, nrObs,  ipoint)); // M x N
		unit<ipoint> MM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef, nrCoef, ipoint)); // M x M
		unit<ipoint> YM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrObs,  1i,     ipoint)); // N x 1
		unit<ipoint> BM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef, 1i,     ipoint)); // M x 1

		attribute<float64> Xt (XtM) := union_data(XtM, const(1.0, Dwelling), Dwelling/area, Dwelling/distance);
		attribute<float64> X  (XM)  := Xt[point_yx(pointcol(id(XM)), pointrow(id(XM)), XtM)];
		attribute<float64> Y  (YM)  := union_data(YM, Dwelling/price);

		attribute<float64> XtX     (MM)   := matr_var(X, MM);
		attribute<float64> XtX_inv (MM)   := matr_inv(XtX);
		attribute<float64> XtY     (BM)   := matr_mul(Xt, Y, BM);
		attribute<float64> Bcol    (BM)   := matr_mul(XtX_inv, XtY, BM);
		attribute<float64> beta    (Coef) := union_data(Coef, Bcol);

		attribute<float64> Yhat      (YM)       := matr_mul(X, Bcol, YM);
		attribute<float64> price_hat (Dwelling) := union_data(Dwelling, Yhat);
		attribute<float64> residual  (Dwelling) := Dwelling/price - price_hat;

		parameter<float64> SSE    := sum(sqr(residual));
		parameter<float64> SST    := sum(sqr(Dwelling/price - mean(Dwelling/price)));
		parameter<float64> R2     := 1.0 - SSE / SST;
		parameter<float64> sigma2 := SSE / float64(nrObs - nrCoef);

		attribute<float64> cov_beta (MM)   := sigma2 * XtX_inv;
		attribute<float64> se       (Coef) := sqrt(cov_beta[point_yx(int32(id(Coef)), int32(id(Coef)), MM)]);
		attribute<float64> t        (Coef) := beta / se;
	}
}

The result:

with $R^2$ = 0.9949 and $s^2$ = 32.0.

Some remarks:

• When $X'X$ is singular, because two regressors are perfectly collinear or there are fewer observations than coefficients, matr_inv results in null in every cell and so does everything computed from it. IsDefined on XtX_inv tells the two cases apart from a data problem.
• Multivariate in the strict sense, several dependent variables at once, needs no new code: stack the dependent variables like the regressors and transpose them into an $N \times K$ matrix Y, make BM an $M \times K$ domain, and $B = (X'X)^{-1}X'Y$ holds every regression in one column.
• Weighted least squares scales the rows of $X$ and $y$ by $\sqrt{w_i}$ through the row relation, as the logistic regression below does.
• For millions of observations X and Xt cost $8 N M$ bytes each. The cells of $X'X$ and $X'y$ are also plain aggregations, sum(x_p * x_q), so with a handful of regressors they can be configured without forming the matrix at all; the matrix route keeps the configuration generic in $M$.

A binary outcome $y_i \in {0, 1}$ with $P_i = P(y_i = 1) = 1 / (1 + e^{-x_i \beta})$ has the log-likelihood $L(\beta) = \sum_i y_i \log P_i + (1 - y_i) \log (1 - P_i)$, with gradient $X'(y - p)$ and Hessian $-X'WX$, where $W$ is the diagonal matrix of $p_i (1 - p_i)$. Newton-Raphson iterates

$$\beta_{k+1} = \beta_k + (X'W_kX)^{-1} X'(y - p_k)$$

and $X'WX = (\sqrt{W}X)'(\sqrt{W}X)$ is matr_var of $X$ with every row scaled by $\sqrt{p_i(1-p_i)}$. One step is a template whose first subitem is the current $\beta$ and whose nextValue is the updated $\beta$; iterate chains the instantiations, feeding each nextValue into the next currValue. From $\beta = 0$, where every $P_i$ is 0.5, the steps converge to machine precision in five to eight iterations for a well-behaved problem.

The example explains whether twelve households moved house from their income and their distance to the city centre.

container regression
{
	unit<uint32> Household : nrofrows = 12
	{
		attribute<float64> income   : [ 22, 28, 31, 35, 40, 44, 47, 52, 58, 63, 70, 78 ];  // k euro per year
		attribute<float64> distance : [  4, 12, 25,  6, 18, 30,  9, 22, 15, 35, 11, 27 ];  // km to the city centre
		attribute<bool>    moved    : [ false, false, false, true, false, true, true, true, true, false, true, true ];
	}

	container binomial_logit
	{
		unit<uint32> Coef : nrofrows = 3 { attribute<string> name : ['const', 'income', 'distance']; }

		parameter<int32> nrObs  := int32(#Household);
		parameter<int32> nrCoef := int32(#Coef);

		unit<ipoint> XM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrObs,  nrCoef, ipoint)) // N x M
		{
			attribute<Household> obs_rel := value(pointrow(id(.)), Household);
		}
		unit<ipoint> XtM := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef, nrObs,  ipoint)); // M x N
		unit<ipoint> MM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef, nrCoef, ipoint)); // M x M
		unit<ipoint> YM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrObs,  1i,     ipoint)); // N x 1
		unit<ipoint> BM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef, 1i,     ipoint)); // M x 1

		attribute<float64> y (Household) := float64(Household/moved);

		attribute<float64> Xt (XtM) := union_data(XtM, const(1.0, Household), Household/income / 10.0, Household/distance / 10.0);
		attribute<float64> X  (XM)  := Xt[point_yx(pointcol(id(XM)), pointrow(id(XM)), XtM)];

		template newton_step
		{
			attribute<float64> currValue (Coef);   // beta at the start of this step

			attribute<float64> Bcol  (BM)        := union_data(BM, currValue);
			attribute<float64> v     (Household) := union_data(Household, matr_mul(X, Bcol, YM));
			attribute<float64> p     (Household) := 1.0 / (1.0 + exp(-v));
			attribute<float64> w     (Household) := p * (1.0 - p);

			attribute<float64> Xw    (XM)   := X * sqrt(w[XM/obs_rel]);
			attribute<float64> H     (MM)   := matr_var(Xw, MM);                            // X'WX
			attribute<float64> H_inv (MM)   := matr_inv(H);
			attribute<float64> G     (BM)   := matr_mul(Xt, union_data(YM, y - p), BM);     // X'(y - p)
			attribute<float64> step  (Coef) := union_data(Coef, matr_mul(H_inv, G, BM));

			attribute<float64> nextValue (Coef) := currValue + step;

			parameter<float64> loglik   := sum(y * log(p) + (1.0 - y) * log(1.0 - p));
			parameter<float64> max_step := max(abs(step));
		}

		unit<uint32> iter : nrofrows = 8 { attribute<string> name := 'step' + string(id(.)); }
		parameter<string> init := 'const(0.0, Coef)';
		container newton := iterate(iter/name, newton_step, init);

		attribute<float64> beta (Coef) := newton/step7/nextValue;
		container final := newton_step(beta);

		attribute<float64> se (Coef) := sqrt(final/H_inv[point_yx(int32(id(Coef)), int32(id(Coef)), MM)]);
		attribute<float64> t  (Coef) := beta / se;

		parameter<float64> loglik  := final/loglik;
		parameter<float64> y_mean  := mean(y);
		parameter<float64> loglik0 := sum(y) * log(y_mean) + (float64(nrObs) - sum(y)) * log(1.0 - y_mean);
		parameter<float64> rho2    := 1.0 - loglik / loglik0;
	}
}

The regressors are scaled to tens of thousands of euros and tens of kilometres, so that the coefficients and the cells of $X'WX$ stay of order one. The result:

The log-likelihood is -8.318 at $\beta = 0$, -5.387 after the first step, -4.879 after the third and -4.877 from the fifth step on; the last step changes no coefficient by more than $10^{-14}$. Against the constant-only log-likelihood of -8.150 this gives McFadden's $\rho^2$ = 0.402.

How the pieces fit:

• iterate copies the template once per name in iter/name and gives the first subitem of copy $k$ the expression step<k-1>/nextValue; the first copy gets the expression in init. The estimate is the nextValue of the last copy, newton/step7/nextValue. The result container also has an item lastValue, but it is untyped, and a typed attribute cannot take it as its calculation rule.
• final := newton_step(beta) is one more case instantiation of the template, at the estimate: it provides $H^{-1}$ for the standard errors, the fitted probabilities final/p and the log-likelihood at the optimum. Its own nextValue is a ninth step that nobody asks for, so it is never calculated.
• Convergence is read from max_step of the last step; when it is not small, enlarge iter. When the log-likelihood keeps rising while a coefficient runs off, the data are (quasi-)completely separated and the maximum likelihood estimate does not exist; drop or merge the regressor that separates the outcomes. Unscaled regressors and poor start values can also make a Newton step overshoot; scale them as above, or use nextValue := currValue + 0.5 * step for the first steps.
• Grouped data, with $n_i$ trials and $y_i$ successes per row, use the same template with w := n * p * (1.0 - p) and y - n * p in the gradient.

Actors $i$ choose from alternatives $j$ with utility $v_{ij} = \sum_p \beta_p X_{ijp}$ and probability $P_{ij} = e^{v_{ij}} / \sum_{j'} e^{v_{ij'}}$; the observed choices $Y_{ij}$ sum to $N_i$ per actor, one for a single recorded choice. Logit regression derives the log-likelihood $L = \sum_{ij} Y_{ij} v_{ij} - \sum_i N_i \log \sum_j e^{v_{ij}}$ and its gradient

$$g_p = \sum_{ij} (Y_{ij} - N_i P_{ij}) X_{ijp}$$

The negative Hessian is a probability-weighted covariance of the regressors within each choice set:

$$-H_{pq} = \sum_i N_i \sum_j P_{ij} (X_{ijp} - \bar X_{ip})(X_{ijq} - \bar X_{iq}), \qquad \bar X_{ip} = \sum_j P_{ij} X_{ijp}$$

so with $Z_{ij,p} = \sqrt{N_i P_{ij}} (X_{ijp} - \bar X_{ip})$ it is $Z'Z$, once more a matr_var, and the Newton step is $(Z'Z)^{-1} g$.

The rows of $X$ are now the pairs of an actor and an alternative: a combine of the actor domain and the alternative domain, with first_rel and second_rel as the relations back. Generic variables such as cost and travel time vary per pair; an alternative-specific constant is a column that is 1 for the pairs of that alternative (form 5 of the specification table on the Logit regression page), and an alternative-specific coefficient is a variable multiplied by such a column (form 3). $\bar X_{ip}$ is a sum of $P_{ij} X_{ijp}$ over the pairs of actor $i$, per coefficient: an $A \times M$ matrix over a domain AM, aggregated through a relation from the cells of XM to the cells of AM.

The example has fifteen travellers choosing between car, bike and transit, with the cost and the travel time of each mode and constants for bike and transit.

container regression
{
	unit<uint32> Traveller : nrofrows = 15;
	unit<uint32> Mode      : nrofrows = 3 { attribute<string> name : ['car', 'bike', 'transit']; }

	unit<uint32> Trip := combine(Traveller, Mode)   // one row per traveller and mode; first_rel: Traveller, second_rel: Mode
	{
		attribute<float64> cost   : [ 5.5, 0, 2.1, 6.4, 0, 3.8, 3.8, 0, 3.9, 8.3, 0, 2, 7.9, 0, 4, 5.3, 0, 2, 4.9, 0, 1.9, 4.9, 0, 1.9, 6.1, 0, 3.4, 7.2, 0, 2.4, 4.3, 0, 1.9, 7.5, 0, 2.8, 3.6, 0, 2.4, 4.7, 0, 3.5, 6.8, 0, 1.8 ];
		attribute<float64> time   : [ 39, 40, 33, 29, 30, 52, 32, 41, 22, 13, 24, 53, 35, 49, 44, 22, 60, 38, 16, 53, 27, 10, 54, 40, 24, 42, 41, 27, 56, 36, 33, 34, 23, 14, 34, 27, 21, 27, 49, 34, 50, 21, 12, 49, 50 ];
		attribute<bool>    chosen : [ false, true, false, false, true, false, false, true, false, false, true, false, true, false, false, false, false, true, false, false, true, true, false, false, false, true, false, false, false, true, false, true, false, false, false, true, true, false, false, false, true, false, true, false, false ];
	}

	container multinomial_logit
	{
		unit<uint32> Coef : nrofrows = 4 { attribute<string> name : ['cost', 'time', 'asc_bike', 'asc_transit']; }

		parameter<int32> nrPairs := int32(#Trip);
		parameter<int32> nrCoef  := int32(#Coef);
		parameter<int32> nrTrav  := int32(#Traveller);

		unit<ipoint> XM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrPairs, nrCoef, ipoint)) // pairs x coefficients
		{
			attribute<Trip>      ij_rel := value(pointrow(id(.)), Trip);
			attribute<Traveller> i_rel  := Trip/first_rel[ij_rel];
			attribute<AM>        am_rel := point_yx(int32(i_rel), pointcol(id(.)), AM);
		}
		unit<ipoint> XtM := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef,  nrPairs, ipoint)); // coefficients x pairs
		unit<ipoint> MM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef,  nrCoef,  ipoint));
		unit<ipoint> RM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrPairs, 1i,      ipoint)); // pairs x 1
		unit<ipoint> BM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrCoef,  1i,      ipoint));
		unit<ipoint> AM  := range(ipoint, point_yx(0i, 0i, ipoint), point_yx(nrTrav,  nrCoef,  ipoint)); // travellers x coefficients

		attribute<float64> y         (Trip)      := float64(Trip/chosen);
		attribute<float64> nrChoices (Traveller) := sum(y, Trip/first_rel);

		attribute<float64> Xt (XtM) := union_data(XtM
			, Trip/cost
			, Trip/time / 10.0
			, float64(Mode/name[Trip/second_rel] == 'bike')
			, float64(Mode/name[Trip/second_rel] == 'transit')
		);
		attribute<float64> X (XM) := Xt[point_yx(pointcol(id(XM)), pointrow(id(XM)), XtM)];

		template newton_step
		{
			attribute<float64> currValue (Coef);

			attribute<float64> Bcol     (BM)        := union_data(BM, currValue);
			attribute<float64> v        (Trip)      := union_data(Trip, matr_mul(X, Bcol, RM));
			attribute<float64> w        (Trip)      := exp(v);
			attribute<float64> denom    (Traveller) := sum(w, Trip/first_rel);
			attribute<float64> p        (Trip)      := w / denom[Trip/first_rel];
			attribute<float64> expected (Trip)      := nrChoices[Trip/first_rel] * p;

			attribute<float64> PX       (XM)   := p[XM/ij_rel] * X;
			attribute<float64> Xbar     (AM)   := sum(PX, XM/am_rel);
			attribute<float64> Z        (XM)   := sqrt(expected[XM/ij_rel]) * (X - Xbar[XM/am_rel]);
			attribute<float64> negH     (MM)   := matr_var(Z, MM);
			attribute<float64> negH_inv (MM)   := matr_inv(negH);
			attribute<float64> G        (BM)   := matr_mul(Xt, union_data(RM, y - expected), BM);
			attribute<float64> step     (Coef) := union_data(Coef, matr_mul(negH_inv, G, BM));

			attribute<float64> nextValue (Coef) := currValue + step;

			parameter<float64> loglik   := sum(y * log(p));
			parameter<float64> max_step := max(abs(step));
		}

		unit<uint32> iter : nrofrows = 8 { attribute<string> name := 'step' + string(id(.)); }
		parameter<string> init := 'const(0.0, Coef)';
		container newton := iterate(iter/name, newton_step, init);

		attribute<float64> beta (Coef) := newton/step7/nextValue;
		container final := newton_step(beta);

		attribute<float64> se (Coef) := sqrt(final/negH_inv[point_yx(int32(id(Coef)), int32(id(Coef)), MM)]);
		attribute<float64> t  (Coef) := beta / se;

		parameter<float64> loglik  := final/loglik;
		parameter<float64> loglik0 := sum(nrChoices) * log(1.0 / float64(#Mode));
		parameter<float64> rho2    := 1.0 - loglik / loglik0;

		attribute<float64> share     (Mode)      := sum(final/p, Trip/second_rel) / float64(nrTrav);
		attribute<Mode>    predicted (Traveller) := Trip/second_rel[max_index(final/p, Trip/first_rel)];
	}
}

The result, converged after six steps:

The log-likelihood is -12.563 against -16.479 for equal shares, $\rho^2$ = 0.238. The predicted shares are 0.267, 0.467 and 0.267, exactly the observed 4, 7 and 4 out of 15: with a constant per alternative the first-order conditions make the predicted counts equal the observed ones. predicted is the most likely mode per traveller. Fifteen observations do not identify two constants well, which the standard errors show; the example is for the mechanics, not for inference.

Variations that need no change to the template:

• Choice sets that differ per actor. Build the pair domain from the available pairs only, for instance with select_with_org_rel on an availability condition or from a spatial join of actors to the alternatives within reach. Every sum is grouped by first_rel, so a shorter choice set simply gets a smaller denominator.
• Sampling of alternatives. When the alternative set is large, draw a sample per actor and add the correction term $\log (1 / q_{ij})$, with $q_{ij}$ the sampling probability, as an offset to v. The estimates stay consistent (McFadden, 1978).
• Aggregate observations. $Y_{ij}$ may be counts, with $N_i$ choices per actor; nrChoices and expected already carry them.
• Alternative-specific coefficients are further columns of Xt: the variable multiplied by the dummy of the alternative.

The nested logit groups the alternatives into nests $k$ and lets the unobserved utilities within a nest correlate. The probability of alternative $j$ in nest $k$ is $P_{ij} = P(j \mid k) , P(k)$ with

$$P(j \mid k) = \frac{e^{v_{ij} / \lambda_k}}{\sum_{j' \in k} e^{v_{ij'} / \lambda_k}}, \qquad IV_{ik} = \log \sum_{j \in k} e^{v_{ij} / \lambda_k}, \qquad P(k) = \frac{e^{\lambda_k IV_{ik}}}{\sum_{k'} e^{\lambda_{k'} IV_{ik'}}}$$

$IV_{ik}$ is the inclusive value, or logsum, of the nest, $\lambda_k \in (0, 1]$ its nesting parameter, and variables of the nest itself can be added to the exponent of $P(k)$. With every $\lambda_k = 1$ the model is the multinomial logit again.

In a configuration the nests are a domain with a relation from the alternatives, and the pairs of an actor and a nest are a second combine. The following fragment computes the probabilities on top of the example above, taking the utilities from the multinomial estimate and putting car and transit in one nest:

container nested_logit
{
	unit<uint32> Nest : nrofrows = 2 { attribute<string> name : ['motorised', 'active']; }
	attribute<Nest> nest_rel (Mode) : [ 0, 1, 0 ];   // car and transit share a nest, bike is alone in its nest

	unit<uint32> TravellerNest := combine(Traveller, Nest);   // first_rel: Traveller, second_rel: Nest
	attribute<TravellerNest> ik_rel (Trip) := combine_data(TravellerNest, Trip/first_rel, nest_rel[Trip/second_rel]);

	attribute<float64> lambda (Nest) : [ 0.7, 1.0 ];   // 1.0 for every nest gives the multinomial logit back

	attribute<float64> v (Trip) := multinomial_logit/final/v;   // the utilities, here from the multinomial estimate

	attribute<float64> w        (Trip)          := exp(v / lambda[nest_rel[Trip/second_rel]]);
	attribute<float64> sum_w    (TravellerNest) := sum(w, ik_rel);
	attribute<float64> IV       (TravellerNest) := log(sum_w);                  // inclusive value of the nest
	attribute<float64> p_within (Trip)          := w / sum_w[ik_rel];           // P(j | nest)

	attribute<float64> wn       (TravellerNest) := exp(lambda[TravellerNest/second_rel] * IV);
	attribute<float64> sum_wn   (Traveller)     := sum(wn, TravellerNest/first_rel);
	attribute<float64> p_nest   (TravellerNest) := wn / sum_wn[TravellerNest/first_rel];   // P(nest)

	attribute<float64> p (Trip) := p_within * p_nest[ik_rel];                   // P(j) = P(j | nest) P(nest)
}

The probabilities sum to one per traveller, and with both $\lambda$ set to 1 they equal multinomial_logit/final/p. Estimating the $\lambda$'s and $\beta$'s together can be done in two ways:

1. Sequentially. The lower level is a multinomial logit of the choice within the chosen nest: the Newton template with its sums grouped by ik_rel instead of first_rel, run on the pairs whose nest was chosen; it estimates $\beta / \lambda$. Then compute $IV_{ik}$ from these estimates, and estimate the upper level as a multinomial logit over the nests with $IV$ and the nest variables as regressors; the coefficient of $IV$ is $\lambda$, a coefficient per nest is form 3 of the specification table. The template serves both levels when it takes the grouping relation as a case parameter. The procedure is consistent but not efficient, and the upper-level standard errors ignore that $IV$ was itself estimated.
2. Full information maximum likelihood. Maximise $L(\beta, \lambda) = \sum_{ij} Y_{ij} \log P_{ij}$ in the same iterate loop. The Hessian is no longer a single weighted cross product, so replace it by the BHHH approximation, the sum over the actors of the outer product of their scores $s_i = \partial L_i / \partial \theta$: with the scores as an $A \times (M + K)$ matrix S, the approximation is matr_var(S, MM) and the step is $(S'S)^{-1} g$. The scores of the nested logit are known in closed form, and can also be approximated by finite differences from instantiations of the probability fragment at $\theta \pm h e_p$. Halve the step when the log-likelihood decreases, and start from the sequential estimates.

A mixed logit draws the coefficients from a distribution and averages the multinomial probability over $R$ draws. In a configuration each draw is an instantiation of the probability fragment with coefficients built from rnd_uniform draws, created with for_each over a domain of draws; the probability is the mean over them and the estimation is again BHHH over the actors' scores. A probit needs the multivariate normal integral and has no counterpart among the GeoDMS functions.

For the multinomial logit the elasticity of $P_{ij}$ to a variable of its own alternative is $\beta_p X_{ijp} (1 - P_{ij})$, and to the same variable of another alternative $j'$ it is $-\beta_p X_{ij'p} P_{ij'}$: attributes of the pair domain that follow from final/p in one line each. The predicted shares and the most likely alternative of the example show how an estimated model is applied; for a new situation, put its pairs and variables in a fresh pair domain and instantiate the template once with the estimated beta.

• Matrix functions, matr_mul, matr_var and matr_inv
• Logit regression for the derivation of the logit likelihood and its specification forms
• iterate, Loop and Template for iteration
• combine, combine_data and union_data for the pair domains and the matrix layout
• Iterative proportional fitting, another estimation loop, and Discrete Allocation, where a logit-like choice is replaced by an allocation under constraints