# How to recover linear probability model for binary classification?

In some binary classification problem I assume that the probability for a positive is exactly linearly dependent on the features $P(y=1|x_1,x_2,x_3,\ldots)=\beta_1x_1+\beta_2x_2+\beta_3x_3+\cdots$. Now for the training I only observe the realization 0 or 1 for $y$. Therefore the probability isn't directly observable, but would rather be some local density (which is hard to estimate for high dimensions).

Is there some type of regression that can determine the coefficients $\beta_i$?

I could try normal regression, but it doesn't seem obvious that any loss function will be able to converge to the same $\beta_i$. What type of regularization and loss function would be most appropriate?

• The specification as you have proposed it is not advisable.. I recommend you look into a generalized linear model (GLM), specifically a logit or probit model. – user44764 Jun 28 '14 at 21:44
• I can look at different models, but I still wouldn't be sure if they can reconstruct the underlying coefficients. Are there theoretical results? Which approach would recover the model? – Gerenuk Jun 28 '14 at 21:50
• GLMs are well-studied. I refer you to Murphy's Probabilistic Machine Learning Chapter 8 for a thorough rundown of logistic regression. – user44764 Jun 28 '14 at 22:03
• What do you mean "converge to the same $\beta$"? As in, converge asymptotically as your sample size grows? Also, yes, a binary/Bernoulli GLM is the go-to parametric approach in this case. A GLM is a consistent estimator for the conditional mean of $y$, which in a Bernoulli distribution is the same as $P(y=1 | \vec{x} )$, so if you use the identity link that's literally what you asked for. You can estimate it with OLS ("regular regression") or maximum likelihood if you want but it's the same model. – shadowtalker Jun 29 '14 at 2:17
• So least squares or just any loss function will always give me the correct model? Certainly one of these loss functions converges faster? – Gerenuk Jun 29 '14 at 16:09

This is the standard Linear Probability Model in the context of Binary Response models. To obtain it, one assumes that there exist an underlying latent/unobservable variable $y^*$, which is linearly dependent on the regressors,

$$y^*=\beta_1x_1+\beta_2x_2+\beta_3x_3+\cdots + u$$

and whose indicator function is the observed $y$, indicating that the latent variable exceeds (or falls below) some threshold (usually zero). In order to arrive at the linear specification for the probability, we have to assume that the error term $u$ in the underlying regression is uniformly distributed, conditional on the regressors.

See this answer for the exact derivation. The model can of course be estimated by least-squares, but has various drawbacks, so frequently a "probit" or "logit" model is used (the first assuming that the underlying error is normally distributed, the second assuming that the underlying error term follows the logistic distribution). These models are usually estimated by maximum likelihood.

After some discussion with the OP and clarification in the comments, it is clear that the OP's framework is represented by the linear probability model where the "error term" is the random draws the OP performs from a $U(0,1)$. Specifically we can postulate the unobservable regression $$y^*_i = -\sum_{j=1}^k\beta_jx_{ji} + u_i,\qquad u_i \sim U(0,1)$$

and the OP determines the values of the observable $y$ by

$$y_i = 1,\;\; \text{if} \;\;u_i \leq \sum_{j=1}^k\beta_jx_{ji} \Rightarrow - \sum_{j=1}^k\beta_jx_{ji} +u_i \leq 0 \Rightarrow y^*_i \leq 0$$

So we have

$$P(y_i = 1 \mid \mathbf X) = P\left(- \sum_{j=1}^k\beta_jx_{ji} +u_i\leq 0\right) =P\left(u_i\leq \sum_{j=1}^k\beta_jx_{ji} \right)$$

But this is the expression for the cumulative distribution function of $u_i$, which is $F_U(z) = P\left(u_i\leq z\right) = z$ for a $U(0,1)$.

In this way we arrive at the equation

$$P(y_i = 1 \mid \mathbf X) = \sum_{j=1}^k\beta_jx_{ji}$$ which is the Linear Probability Model, on firm stochastic ground. Note that now the "error term" is no longer "present" -but it has determined the structure of the equation, through its stochastic properties.

Estimating the $\beta$'s from this last equation by Least Squares does not require specifying some additional "error term" -it has now become more of a mathematical approximation endeavor, rather than a stochastic regression. Another way of course would be to estimate the model by maximum likelihood. Here the distributional specification on $y$ is "safe": it is by construction a Bernoulli random variable, whose conditional parameter $p$ is determined by $\sum_{j=1}^k\beta_jx_{ji}$. We can write the Bernouli probability mass function as

$$f_Y(y_i\mid \mathbf X) = p_i^{y_i}(1-p_i)^{1-y_i} = \left(\sum_{j=1}^k\beta_jx_{ji}\right)^{y_i}\left(1-\sum_{j=_1}^k\beta_jx_{ji}\right)^{1-y_i}$$

from which we can proceed to form the log-likelihood of the sample and perform maximum likelihood.

Estimator properties and other related issues to the estimation are easily accessible.

• +1 For giving the answer the OP needs, but I think he's more interested in guarantees/proofs that the coefficients are stable and that some kind of meaningful loss will function will allow him to recover the same coefficients.. or at least that's my interpretation of the question. – user44764 Jun 29 '14 at 2:07
• How come I need a new error term $u$ and some distribution for it? I'm not too familiar with all terminology, but from my understanding I imagine simulating it by generating latent probabilities by a linear model and then using these probabilities to binarize them to 0 or 1. There is no place for an error term. Or equivalently the error term is $P(u=-p)=1-p$, $P(u=1-p)=p$ where $p=\sum \beta_ix_i$. I have the impression that this is different from your model? – Gerenuk Jun 29 '14 at 12:37
• @Gerenuk You start by a conditional probability on the left-hand-side. Probabilities are expressed by distribution functions. You condition on the $x$'s, so the conditional probability/distribution of $Y$ cannot involve the distributions of the $x$'s. But it must come from somewhere. There must be something else that is random, otherwise this could not be a probability. And this something else is the underlying error term. You write "generate latent probabilities". How are you going to generate them without assuming a distribution? And if you assume a distribution, whose distribution is it? – Alecos Papadopoulos Jun 29 '14 at 19:14
• Maybe I'm not using the right statistics terms, so I try to explain what a computer program would do. The generator selects some $\beta_i$ and calculates $p(x)=\sum \beta_ix_i$. For each $x$ it will use $y=1$ with that probability $p(x)$ (I run take a uniform random number $[0;1]$ and use $p(x)$ as a threshold for $y=1$). Now my second model fitting program is expected to take all $y(x)$ as input and yield an good estimates of $\beta_i$ given a limited amount of data. – Gerenuk Jun 30 '14 at 6:19
• You mean that you draw random numbers from a uniform $U(0,1)$ (denote one of these random numbers $u$), and if $u$ is smaller or equal to $\sum\beta_ix_i$ you assign the value $1$ to $y$, and you set $y=0$ otherwise? – Alecos Papadopoulos Jun 30 '14 at 9:08