8
$\begingroup$

I want to fit a curve $f(x) = mx+b$ on my data points $x_1, \ldots, x_N$ using linear regression with a single predictor.

However, the cost function is not even, rather, it has different weights on each side, i.e.:

$$ E = \frac{1}{N}\sum_{i=1}^N{\text{Cost}(f(x_i) - y_i)} . $$

where: $$\text{Cost}(v) = \begin{cases}v^2 &v\leq 0 \\\alpha \cdot v^2 & v> 0\end{cases} $$$$ \alpha > 0 $$

Are there any well-known methods for finding such a line for an arbitrary $\alpha$ value

I specifically wonder how should I find the line $mx+b$ that is totally under the points? (i.e. $\alpha=\infty$)

$\endgroup$
6
  • 7
    $\begingroup$ I believe this is what is called asymmetric least squares or expectiile regression. See, for instance, freakonometrics.hypotheses.org/files/2017/05/erasmus-1.pdf . Perhaps the R package expectreg an handle this, $\endgroup$ Commented Jul 9, 2018 at 0:37
  • $\begingroup$ I am also interested. What packages are there for doing local (polynomial) quantile and expectile regressions in R? $\endgroup$
    – Viktor
    Commented Jul 9, 2018 at 2:02
  • $\begingroup$ Re your last question: simply choose any $(m,b)$ that works! For instance, $m=0$ and $b=\min(y_i)$ obviously qualifies. Concerning the general question: first apply Calculus to solve the problem for the model $f(x)=b$ so you can understand the nuances. $\endgroup$
    – whuber
    Commented Jul 9, 2018 at 12:34
  • 1
    $\begingroup$ @whuber That is indeed not correct, because we still need to minimize Cost(v) on for positive error values. $\endgroup$
    – Ali
    Commented Jul 9, 2018 at 14:00
  • 1
    $\begingroup$ @whuber Thanks. Instead of setting $\alpha=\infty$, can we convert it into a constraint on the function? e.g.: Minimize $\sum{Cost((f(x_i)-y_i)^2)}$ s.t. $f(x_i) \leq y_i$ $\endgroup$
    – Ali
    Commented Jul 10, 2018 at 10:25

2 Answers 2

1
$\begingroup$

As Mark L. Stone comments, your more general problem is indeed known as expectile regression. A nice overview can be found in Waltrup et al. (2015). We need to map your weights of $1$ and $\alpha$ to expectiles that sum to $1$; a little algebra shows that you need $\frac{1}{1+\alpha}$-expectiles.

Again per Mark's comment, the expectreg package for R will do the fitting for you (and even include penalization if you want). Here are some random data with a correlation and the expectile regression fits for $\alpha\in\{0.1, 1, 10\}$:

expectile regression 1

R code:

n_points <- 25

set.seed(1)
xx <- rnorm(n_points)
yy <- xx+rnorm(n_points)
dataset <- data.frame(x=xx,y=yy)
x_predict <- range(xx)

alpha <- c(0.1,1,10)
expectiles <- 1/(1+alpha)

library(expectreg)
opar <- par(mai=c(.8,.8,.1,.1))
    plot(dataset,las=1,pch=19)
    for ( ii in seq_along(alpha) ) {
        model <- expectreg.ls(y~x,dataset,lambda=0,expectiles=expectiles[ii])
        lines(x_predict,predict(model,newdata=data.frame(x=x_predict))$fitted,col=ii,lwd=2)
    }
    legend("bottomright",lwd=2,col=seq_along(alpha),legend=paste("alpha =",alpha))
par(opar)

Now, for your more specific problem of finding a line that passes under all points with minimal sum of squared errors. There are multiple ways to go about this:

  1. In the setting of expectile regression, you can simply work with the $0$-expectile (which is simply the limit, $\lim_{\alpha\to\infty}\frac{1}{1+\alpha}=0$). Just use the parameter expectiles=0 in the call to expectreg.ls() above, and you get your line. Symmetrically, expectiles=1 will give you a line above your points (and a ton of warnings, which don't look very serious to me).

  2. Alternatively, finding a line below (or above) all points with minimal costs is a straightforward quadratic programming exercise. Your objective function is the sum of squared residuals, which you want to minimize. Your constraints are that $ax_i+b\leq y_i$ for all $i$, and these are obviously linear. You can just feed this into the quadratic programming solver of your choice.

  3. Finally, a geometric approach would be to compute the convex hull of your point cloud. A line that is completely below your point cloud is determined by one of the finitely many segments of this hull below the cloud. So you can simply look at each such segment, extend it to a line in both directions, and compute the sum of squared residuals - and finally use the line with the smallest such sum. Again, the same works to find a minimal squared residual line above your point cloud.

$\endgroup$
-3
$\begingroup$

Here is example code using the Python scipy.optimize.differential_evolution genetic algorithm module implementing a "brick wall" that gives a very large error if the genetic algorithm finds parameters that yield any predicted value above or below that of any data point per a code switch. It works in my tests when I flip the "upper/lower" switch in the code.

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import warnings

from scipy.optimize import differential_evolution


xData = np.array([5.0, 6.1, 7.2, 8.3, 9.4])
yData = np.array([ 10.0,  18.4,  20.8,  23.2,  35.0])


def func(data, a, b):
    return a * data + b


# function for genetic algorithm to minimize (sum of squared error)
# this contains the "brick wall" switch for upper/lower
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    for i in range(len(val)):
        if val[i] < yData[i]: # ****** upper/lower switch ******
            val[i] = 1.0E10
    return np.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)

    parameterBounds = []
    maxSlope = (maxY - minY) / (maxX / minX)
    parameterBounds.append([-maxSlope, maxSlope]) # parameter bounds for a
    parameterBounds.append([-maxY, maxY]) # parameter bounds for b

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# generate initial parameter values
geneticParameters = generate_Initial_Parameters()

# create values for display of fitted peak function
a, b = geneticParameters
y_fit = func(xData, a, b)

plt.plot(xData, yData, 'D') # plot the raw data
plt.plot(xData, y_fit) # plot the equation using the fitted parameters
plt.show()

print('parameters:', geneticParameters)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.