Maximum likelihood estimation

Suppose data set ${\displaystyle \{y_{i},x_{i}\}_{i=1}^{n}}$ contains n independent statistical units corresponding to the model above.

For the single observation, conditional on the vector of inputs of that observation, we have:

${\displaystyle P(y_{i}=1|x_{i})=\Phi (x_{i}'\beta )}$^{[clarification needed]}

${\displaystyle P(y_{i}=0|x_{i})=1-\Phi (x_{i}'\beta )}$

where ${\displaystyle x_{i}}$ is a vector of ${\displaystyle K\times 1}$ inputs, and ${\displaystyle \beta }$ is a ${\displaystyle K\times 1}$ vector of coefficients.

The likelihood of a single observation ${\displaystyle (y_{i},x_{i})}$ is then

${\displaystyle {\mathcal {L}}(\beta ;y_{i},x_{i})=\Phi (x_{i}^{\operatorname {T} }\beta )^{y_{i}}[1-\Phi (x_{i}^{\operatorname {T} }\beta )]^{(1-y_{i})}}$

In fact, if ${\displaystyle y_{i}=1}$, then ${\displaystyle {\mathcal {L}}(\beta ;y_{i},x_{i})=\Phi (x_{i}^{\operatorname {T} }\beta )}$, and if ${\displaystyle y_{i}=0}$, then ${\displaystyle {\mathcal {L}}(\beta ;y_{i},x_{i})=1-\Phi (x_{i}^{\operatorname {T} }\beta )}$.

Since the observations are independent and identically distributed, then the likelihood of the entire sample, or the joint likelihood, will be equal to the product of the likelihoods of the single observations:

${\displaystyle {\mathcal {L}}(\beta ;Y,X)=\prod _{i=1}^{n}\left(\Phi (x_{i}^{\operatorname {T} }\beta )^{y_{i}}[1-\Phi (x_{i}^{\operatorname {T} }\beta )]^{(1-y_{i})}\right)}$

The joint log-likelihood function is thus

${\displaystyle \ln {\mathcal {L}}(\beta ;Y,X)=\sum _{i=1}^{n}{\bigg (}y_{i}\ln \Phi (x_{i}^{\operatorname {T} }\beta )+(1-y_{i})\ln \!{\big (}1-\Phi (x_{i}^{\operatorname {T} }\beta ){\big )}{\bigg )}}$

The estimator ${\displaystyle {\hat {\beta }}}$ which maximizes this function will be consistent, asymptotically normal and efficient provided that ${\displaystyle \operatorname {E} [XX^{\operatorname {T} }]}$ exists and is not singular. It can be shown that this log-likelihood function is globally concave in ${\displaystyle \beta }$, and therefore standard numerical algorithms for optimization will converge rapidly to the unique maximum.

Asymptotic distribution for ${\displaystyle {\hat {\beta }}}$ is given by

${\displaystyle {\sqrt {n}}({\hat {\beta }}-\beta )\ {\xrightarrow {d}}\ {\mathcal {N}}(0,\,\Omega ^{-1}),}$

where

${\displaystyle \Omega =\operatorname {E} {\bigg [}{\frac {\varphi ^{2}(X^{\operatorname {T} }\beta )}{\Phi (X^{\operatorname {T} }\beta )(1-\Phi (X^{\operatorname {T} }\beta ))}}XX^{\operatorname {T} }{\bigg ]},\qquad {\hat {\Omega }}={\frac {1}{n}}\sum _{i=1}^{n}{\frac {\varphi ^{2}(x_{i}^{\operatorname {T} }{\hat {\beta }})}{\Phi (x_{i}^{\operatorname {T} }{\hat {\beta }})(1-\Phi (x_{i}^{\operatorname {T} }{\hat {\beta }}))}}x_{i}x_{i}^{\operatorname {T} },}$^{[citation needed]}

and ${\displaystyle \varphi =\Phi '}$ is the Probability Density Function (PDF) of standard normal distribution.

Semi-parametric and non-parametric maximum likelihood methods for probit-type and other related models are also available.^[4]

Berkson's minimum chi-square method

This method can be applied only when there are many observations of response variable ${\displaystyle y_{i}}$ having the same value of the vector of regressors ${\displaystyle x_{i}}$ (such situation may be referred to as "many observations per cell"). More specifically, the model can be formulated as follows.

Suppose among n observations ${\displaystyle \{y_{i},x_{i}\}_{i=1}^{n}}$ there are only T distinct values of the regressors, which can be denoted as ${\displaystyle \{x_{(1)},\ldots ,x_{(T)}\}}$. Let ${\displaystyle n_{t}}$ be the number of observations with ${\displaystyle x_{i}=x_{(t)},}$ and ${\displaystyle r_{t}}$ the number of such observations with ${\displaystyle y_{i}=1}$. We assume that there are indeed "many" observations per each "cell": for each ${\displaystyle t,\lim _{n\rightarrow \infty }n_{t}/n=c_{t}>0}$.

Denote

${\displaystyle {\hat {p}}_{t}=r_{t}/n_{t}}$

${\displaystyle {\hat {\sigma }}_{t}^{2}={\frac {1}{n_{t}}}{\frac {{\hat {p}}_{t}(1-{\hat {p}}_{t})}{\varphi ^{2}{\big (}\Phi ^{-1}({\hat {p}}_{t}){\big )}}}}$

Then Berkson's minimum chi-square estimator is a generalized least squares estimator in a regression of ${\displaystyle \Phi ^{-1}({\hat {p}}_{t})}$ on ${\displaystyle x_{(t)}}$ with weights ${\displaystyle {\hat {\sigma }}_{t}^{-2}}$:

${\displaystyle {\hat {\beta }}={\Bigg (}\sum _{t=1}^{T}{\hat {\sigma }}_{t}^{-2}x_{(t)}x_{(t)}^{\operatorname {T} }{\Bigg )}^{-1}\sum _{t=1}^{T}{\hat {\sigma }}_{t}^{-2}x_{(t)}\Phi ^{-1}({\hat {p}}_{t})}$

It can be shown that this estimator is consistent (as n→∞ and T fixed), asymptotically normal and efficient.^{[citation needed]} Its advantage is the presence of a closed-form formula for the estimator. However, it is only meaningful to carry out this analysis when individual observations are not available, only their aggregated counts ${\displaystyle r_{t}}$, ${\displaystyle n_{t}}$, and ${\displaystyle x_{(t)}}$ (for example in the analysis of voting behavior).

Gibbs sampling

Gibbs sampling of a probit model is possible because regression models typically use normal prior distributions over the weights, and this distribution is conjugate with the normal distribution of the errors (and hence of the latent variables Y^*). The model can be described as

${\displaystyle {\begin{aligned}{\boldsymbol {\beta }}&\sim {\mathcal {N}}(\mathbf {b} _{0},\mathbf {B} _{0})\\[3pt]y_{i}^{\ast }\mid \mathbf {x} _{i},{\boldsymbol {\beta }}&\sim {\mathcal {N}}(\mathbf {x} _{i}^{\operatorname {T} }{\boldsymbol {\beta }},1)\\[3pt]y_{i}&={\begin{cases}1&{\text{if }}y_{i}^{\ast }>0\\0&{\text{otherwise}}\end{cases}}\end{aligned}}}$

From this, we can determine the full conditional densities needed:

${\displaystyle {\begin{aligned}\mathbf {B} &=(\mathbf {B} _{0}^{-1}+\mathbf {X} ^{\operatorname {T} }\mathbf {X} )^{-1}\\[3pt]{\boldsymbol {\beta }}\mid \mathbf {y} ^{\ast }&\sim {\mathcal {N}}(\mathbf {B} (\mathbf {B} _{0}^{-1}\mathbf {b} _{0}+\mathbf {X} ^{\operatorname {T} }\mathbf {y} ^{\ast }),\mathbf {B} )\\[3pt]y_{i}^{\ast }\mid y_{i}=0,\mathbf {x} _{i},{\boldsymbol {\beta }}&\sim {\mathcal {N}}(\mathbf {x} _{i}^{\operatorname {T} }{\boldsymbol {\beta }},1)[y_{i}^{\ast }<0]\\[3pt]y_{i}^{\ast }\mid y_{i}=1,\mathbf {x} _{i},{\boldsymbol {\beta }}&\sim {\mathcal {N}}(\mathbf {x} _{i}^{\operatorname {T} }{\boldsymbol {\beta }},1)[y_{i}^{\ast }\geq 0]\end{aligned}}}$

The result for ${\displaystyle {\boldsymbol {\beta }}}$ is given in the article on Bayesian linear regression, although specified with different notation.

The only trickiness is in the last two equations. The notation ${\displaystyle [y_{i}^{\ast }<0]}$ is the Iverson bracket, sometimes written ${\displaystyle {\mathcal {I}}(y_{i}^{\ast }<0)}$ or similar. It indicates that the distribution must be truncated within the given range, and rescaled appropriately. In this particular case, a truncated normal distribution arises. Sampling from this distribution depends on how much is truncated. If a large fraction of the original mass remains, sampling can be easily done with rejection sampling—simply sample a number from the non-truncated distribution, and reject it if it falls outside the restriction imposed by the truncation. If sampling from only a small fraction of the original mass, however (e.g. if sampling from one of the tails of the normal distribution—for example if ${\displaystyle \mathbf {x} _{i}^{\operatorname {T} }{\boldsymbol {\beta }}}$ is around 3 or more, and a negative sample is desired), then this will be inefficient and it becomes necessary to fall back on other sampling algorithms. General sampling from the truncated normal can be achieved using approximations to the normal CDF and the probit function, and R has a function rtnorm() for generating truncated-normal samples.