# How best to estimate regression parameters subject to constraints?

### The setting

Consider a least squares model for $$y$$ as a function of $$x,$$ possibly nonlinear in the parameters. Abstractly this can be expressed as

$$y = f(x;\theta) + \varepsilon$$

with the usual assumptions: the $$\varepsilon$$ are iid zero-mean errors of finite variance $$\sigma^2.$$

Further suppose this model represents a well-established mathematical or physical law that constrains the "reasonable" values of the parameter $$\theta.$$

As a concrete example, consider a standard model in radiological measurements such as a (simplified) Lorentzian peak at $$x=0$$ added to a constant background:

$$f(x,(\alpha,\beta,\gamma)) = \frac{\alpha}{\beta + x^2} + \gamma.$$

All three parameters $$\theta = (\alpha,\beta,\gamma)$$ are of interest and should be estimated as accurately as possible. However, the only physically realistic values of the parameters are non-negative (and $$\beta$$ is nonzero). Here is an illustration.

Generally, suppose there is a known nontrivial set $$\Theta$$ of realistic values of $$\theta.$$ I will refer to $$\Theta$$ as the constraints.

A fundamental question is whether the peak is present at all. That can be tested with the null hypothesis $$H_0:\alpha=0,$$ a situation that often holds in reality (maybe the peak is produced by a contaminant that usually is not present, for instance).

When fitting such a model, it is possible for any of the parameter estimates $$\hat\theta$$ to lie outside the constraints. (That is, $$\hat\theta\notin \Theta.$$) For instance, $$\hat\alpha$$ might be negative when (a) the true value $$\alpha$$ is much smaller than the error standard deviation $$\sigma$$ and (b) due to random (bad) luck, it looks like the peak is inverted. This happens frequently (about half the time) when $$\alpha$$ truly is zero.

### The question

One has two options in these circumstances:

1. Use whatever estimates $$\hat\theta$$ the least squares procedure provides, even if they lie outside the constraints.

2. Force the least squares solution to conform to the constraints.

There are cogent arguments for and against both options: (1) can yield unrealistic estimates, but when applied to a large number of datasets might provide better average estimates on the whole; whereas (2) is forced to yield realistic estimates, but that can introduce a bias. For instance, forcing $$\alpha$$ to be nonnegative will produce positive nonzero estimates from time to time even when $$\alpha$$ is truly zero. In general, we're considering the possibility $$H_0: \theta\in\partial\Theta$$ that the parameter lies on the boundary of the constraints and wish to test the alternative $$H_A:\theta\in \Theta \setminus \partial\Theta$$ (the interior of the constraint space).

Is one of these options preferable? If so, which one and why? Is one of these options so problematic that it should not be used? If so, which one and why?

I also invite suggestions for additional options (such as regularization, penalized least squares, and so on), provided they can be shown superior to both (1) and (2).

I welcome answers from varying perspectives (classical, Bayes, etc.) and have enough bounty points to reward multiple good answers ;-).

• I would submit that another advantage of approach 1) is that getting unrealistic estimates of model parameters can prompt us to take a step back from the modeling phase of data analysis and head back to the EDA phase, or maybe even the data collection phase, because it allows our expectations to be contradicted if something unexpected is going on. Commented Oct 11, 2023 at 17:40
• My concern about option #1 is what to do with a preposterous value. For instance, we could do something wild and get a variance estimate below zero. Okay, we have our $\hat\sigma^2_{\text{weird}}<0$, but what does that mean?
– Dave
Commented Oct 11, 2023 at 17:42
• In empirical likelihood estimation it is all about incorporating constraints en.m.wikipedia.org/wiki/Empirical_likelihood maybe a potential answer can discuss this method as well, although this may not be needed if the distribution of measurements is known. Incorporating constraints may be required when convergence to a feasible mimimum of the loss function is a problem. Commented Oct 11, 2023 at 17:53
• @John Agreed, but this is not EDA or a modeling phase. Imagine that this fitting procedure is, say, converting an electrical response in a consumer scale to a weight and it will be used in an automatic turnkey capacity to drive an LED readout of that weight.
– whuber
Commented Oct 11, 2023 at 18:39
• This seems to be a very cool example of the bias-variance trade-off Hastie, et al discuss with regularization methods. My question to you is whether you’re more interested in accurate information or with consistent information. Obviously both solutions you pose are feasible, it sounds like you’re interested in a more philosophical answer, yes? Commented Oct 12, 2023 at 17:17

Is one of these options preferable? If so, which one and why? Is one of these options so problematic that it should not be used? If so, which one and why?

Depends on your goal. If you trust your priors, model specification and error distribution, you would like to include as much of those prior beliefs in the model as you can.

If you want to perform an analysis on a large number of samples, on the other hand, you might want to perform a different analysis altogether. Of course, if you need to find out whether there is some systematic bias in your measurements, you need to use an estimator that has minimal bias -- but then could choose to model the bias directly.

So in general, I would dismiss your pro towards the method (1): can yield unrealistic estimates, but when applied to a large number of datasets might provide better average estimates on the whole, and only worry about the bias if it is the quantity you want to model.

## Frequentist vs. Bayesian interpretation

There was a paper I found comparing frequentist constraints and bayesian priors, the difference of which comes down to interpretation as well as practical usage. Since the frequentist approach treats the parameter $$\theta$$ as unknown, but fixed, element of $$\Theta$$, it does not make much sense to me to allow values that are physically impossible -- i.e. there is no sensible interpretation for those. In practice, then, this means that the estimation must be restricted to a set of parameters that satisfy the constraints, if you wish to set them. Alternatively, you simply accept that your model doesn't actually model the phenomena.

In the Bayesian framework, the parameter itself is a random variable, for which we know at least that $$\theta \in \Theta$$. Of course, there are infinitely many distributions that can be used to assign total probability of 1 to that space. The interpretation of prior distributions to restrict parameter values is straightforward, and you get the (posterior) probability distribution of the parameter values, which, again, has quite simple interpretation. It could include physically impossible values, but in that case it is clear that they were, also, included in your prior beliefs.

I would say, then, that in the Bayesian case it is quite OK to not restrict the parameter to a certain range, and this has a simple explanation: We chose to allow the parameter to take those values, even if we knew that they were impossible. Conversely if you choose to restrict the parameter space by using a prior, the meaning of that is clear and reflects the analyst's prior beliefs: There is zero probability that the parameter could take these values.

If you are concerned about systematic biases introduced by the prior, one option would be to use cross-validation to check the distribution of the errors.

## Bayesian approach in practice

Most common (and trivial) case where you restrict the parameter space by a prior functions is when setting a prior for variance of normal distribution: Inverse gamma distribution is often used. For most parameters there is often enough data that accurate prior distributions matter little, and noninformative, computationally easy (e.g. conjugate) priors are often chosen.

More relevant to your question, perhaps, is this is done with physical models in practice. One paper that might be of interest is A Bayesian Search for the Higgs Particle, where one distribution of the prior is based on the theoretical, expected mass of the Higgs Boson -- and another on what we might see if it does not exist.

The same author also had another paper that is more general and thus maybe more relevant to your question.

Option 1 is better. Here's why:

My favorite professor (Herman Friedman) used to say if you're not surprised, you haven't learned anything.''

Option 2 is a way of telling the data not to surprise you. It's the antithesis of science.

But what if option 1 gives you surprising values? Then you figure out why. You learn something. Maybe it's something mundane (often, data entry errors or coding errors) and you just fix the problem. But once in a wonderful while, it isn't mundane and you discover something. The model is wrong. Or, at least, wrong in your application.

Suppose you are interested in the heights of men and women. Maybe you have some theory about covariates. But you know men are taller than women, so you constrain the estimates. But what if you discovered some subgroup of people where men are not taller? That's really, really interesting. And someone else will find it, not you, and there goes your great discovery.

(I don't understand your example well enough to use it, but if it is really an absolute physical constraint, then you are in the first paragraph, and finding the weird estimates lets you discover that.)

• Isn’t your argument against option #2 equivalent to arguing against informative prior distributions?
– Dave
Commented Oct 11, 2023 at 18:30
• In the circumstances you posit, I agree, But I guess I need to re-emphasize that having a negative peak is not a surprise: it's non-physical. This isn't a matter of scientific discovery, but purely of estimation and testing. Here's another example: suppose the model estimates weights of solid objects placed on a scale. Are you arguing that the scale should report negative weights to its users?
– whuber
Commented Oct 11, 2023 at 18:37
• @Dave No. You can have prior beliefs, whether you formalize them in a Bayesian framework or not. But you can still be surprised. Bayesians let the informative priors be adjusted by data. Commented Oct 11, 2023 at 18:49
• I understand where you're coming from, and am sympathetic, but I have to disagree with your logic. First, there's a question of how significant these impossible results are: they might be within a couple standard errors of the parameter boundary. That's no evidence of bad data or bad code. Second, if they are way off, then indeed that's a nice check of the algorithm, but it doesn't really answer the question. What should one do when the code and the data are good?
– whuber
Commented Oct 11, 2023 at 18:57
• Good point. My main reason for thinking Bayesian for this problem is that it doesn't require a dichotomous solution (hard constraints). Favoritism of the constraint still allows some possibility of negative values, for example. Commented Oct 16, 2023 at 19:23

During reading What are some motivations for using nonnegative least squares? I came across nonnegative least squares which imposes the constraint of all coefficients being nonnegative. The main argument seems to be that prior knowledge is encoded in the constraints and solutions which violate the constraints are known to be unreasonable.

• Thank you. This restates point (2) in my question -- but it doesn't seem to take us any closer to a resolution of the tension between points (1) and (2).
– whuber
Commented Oct 14, 2023 at 14:54

# One dimensional simple model

A simple model that illustrates a bit what is going one is a single parameter model

$$X_i \sim N(\theta,1) \quad \text{with the constraint \theta \geq 0}$$

For a sample of size $$n$$

• the sufficient statistic is $$T=\bar{X}$$ (and we could use this as an unconstrained estimate $$\hat{\theta}_{unconstrained} = T$$)
• the constrained maximum likelihood estimate is $$\hat{\theta}_{constrained} = \max(0,\bar{X})$$.

Histograms for the distributions for these, when $$n=5$$ and $$\theta = 0.5$$, look like below. The distribution of the constrained estimate is a rectified Gaussian distribution.

That peak at $$mle=0$$ is a point mass and represents the tail from $$T<0$$ that is now closer to the true value $$\theta$$. In this case the boundaries are making the mean squared error smaller.

### Variance, bias, and combining multiple experiments

The sufficient statistic $$T$$ has no bias but the variance is large and the mean squared error will be larger.

$$\text{mean}(T) = \theta\\ \text{var}(T) = \frac{1}{n}\\ MSE(T) = \frac{1}{n}$$

The constrained maximum likelihood estimate $$\hat{\theta}_{constrained}$$ has a bias, but due to a smaller variance the mean squared error will be smaller.

(for derivations see: Expectation and Variance of Gaussian going through Rectified Linear or Sigmoid function)

$$\text{mean}(\hat{\theta}_{constrained}) = \theta (1-p) + \frac{\phi(z)}{\sqrt{n}}\\ \text{var}(\hat{\theta}_{constrained}) = \theta^2 p(1-p) + \frac{\theta}{\sqrt{n}} \phi(z)(2p-1) + \frac{1}{n} (1-p-\phi(z)^2) \\ MSE(\hat{\theta}_{constrained}) = \frac{1}{n} (1-p) - \frac{\theta}{\sqrt{n}} \phi(z) + \theta^2 p$$

Where $$z = -\theta \sqrt{n}$$ and $$p = \Phi(z)$$ and $$\phi$$ and $$\Phi$$ are the standard Gaussian pdf and cdf.

The issue/consideration for the one dimensional simple model is

(1) can yield unrealistic estimates, but when applied to a large number of datasets might provide better average estimates on the whole;

When we have multiple estimates then the variance will reduce and the error of the average of estimates will approach the bias.

So when we use averages of estimates then using biased estimates can be a disadvantage, albeit these estimates having a smaller error when compared individually.

Depending on how the model is used we may prefer methods 1 or 2.

• If we are not using averages then method 2 is better because the mean squared error is larger.
• If we will use some averaging of several estimates (and we are not able to use all the data at once making a single estimate instead of averaging multiple estimates) then model 1 might be better. It will depend on the ratio between variance and bias.

Method 2 might give a better estimate, but that estimate might not be the best statistic for further analysis. In the example above this has been made more salient by method 1 being equivalent to the sufficient statistic, and method 2 is effectively loosing information. The same principle might be true for other methods (like regularised regression) where some estimate has 'improved' properties, but results in a loss of information.

• At the end I wrote "It will depend on the ratio between variance and bias". In a situation like non-negative least squares with many correlated components as in the elsewhere linked question What are some motivations for using nonnegative least squares?, the variance of the unconstrained model will be very large. In such a case, even when method 1 may lead to a better estimate for averages, it only does so when the number of averaged samples is extremely large. Commented Jan 18 at 10:50
• In addition, in the case of peak detection and the prior assumption that peaks are mostly absent, it might not be useful to allow negative values. These negative values will coincide with other values being positive. An interesting alternative for nnls might be LASSO stats.stackexchange.com/a/426984 Commented Jan 18 at 10:52
• Thank you for your analysis, +1. But please bear in mind the "fundamental question," that is, whether the peak is present at all. That is so very important (it's the question of whether a substance has been detected and that can have huge consequences) that variance isn't even a secondary consideration except insofar as it might be a proxy for estimation precision.
– whuber
Commented Jan 18 at 14:47
• @whuber the example here is a lot more simple to investigate what is going on. I believe that this main issue is that approach 2 might be a good estimate, but it is not a good (informative) statistic. If applied to some hypothesis test, then a similar issue should be going on, that we rather have as data the different values T instead of the different values $\hat{\theta}_{constrained}$. Commented Jan 18 at 15:12