Clipped univariate OLS regression I have $x, y$ data, and would like to compute an OLS best fit curve of the form:
$$\hat{y} = m*\mathrm{clip}(x, x_{min}, x_{max}) + b$$
Here, $\mathrm{clip}()$ refers to the numpy clip() function, and $m, b, x_{min}, x_{max}$ are the parameters to optimize.
Is there a package (pyearth?) that can solve this problem? Alternatively, is there an efficient algorithm one can implement to solve this problem? I am ok with exact or approximate solutions.
What if we make it multivariate instead of univariate?
 A: If you're open to using R, this would be the mcp model:
model = list(
  y ~ 1,  # b
  ~ 0 + m,  # b + add slope
  ~ 0 # flat line from here
)
fit = mcp(model, data)

mcp will estimate the common intercept $b$, the slope $m$, and the two change points ($x_{min}$, $x_{max}$) for you. Read more at the mcp website.
As of version 0.3, mcp doesn't support multivariate data, unfortunately.
A: What you ask for is an optimization algorithm and although I don't know about the theoretical guarantees for your case you could use a [Gradient Descent algorithm][1]. This is not really an ad-hoc solution but it is a good start, you can see if it works. To compute this, you only need to compute the gradient, the clip part is not differentiable but let us take the almost everywhere derivative and try with that. For simplicity I will only treat the case where $x_{min}=x_{max}=c$, the general case can be inferred from here.
Let $\psi_c(x)=x 1\{ |x|\le c\}+c\mathrm{sign}(x) 1\{|x|>c\}$ where $1$ is the indicator function (in higher dimension this could be replaced with $\psi_c^{multi}(x)=\frac{x}{\|x\|}\psi_c(\|x\|)$ for example). If I am not mistaken, your problem is to find
$$(a,b,c) \in \arg\min J(a,b,c)$$
where
$$ J(a,b,c)=\frac{1}{n}\sum_{i=1}^n \left(y_i-(a\psi_c(x_i)+b)\right)^2 $$
Then, the derivatives of $J$ with respect to $a$ and $b$ are the same as in OLS, and the (almost everywhere) derivative with respect to $c$ is
$$\frac{d}{dc}J(a,b,c)=-\frac{2}{n}\sum_{i=1}^n \left(y_i-(a\psi_c(x_i)+b)\right)a\mathrm{sign}(x_i)1\{c \le |x_i|\}. $$
and then, you can update $c$ using this gradient as you would do with an ordinary gradient descent except that you the maximum with zero if it happens to be negative at a time (this is a projected gradient step). In practice, this gradient is a correction on the points that are clipped by your function.
EDIT : Here is a code. It is very preliminary, my gradient descent is not well tuned and I advise you to use a line search or to use more powerful algorithm.
import numpy as np
import matplotlib.pyplot as plt

def psi(x,c):
    return x*(np.abs(x)<c)+c*(2*(x>0)-1)*(np.abs(x)>c)

X = np.random.standard_t(3, size=100) # Heavy-tailed with variance 3
y =  5*psi(X,0.5)+2

# Initialization
a = 1
b = 1
c = 1
M = 1000 # number of steps
step_size = 0.05

# Gradient descent
for f in range(M):
    grad = [-2*np.mean(psi(X,c)*(y-a*psi(X,c)-b)), -2*np.mean(y-a*psi(X,c)-b),
            -2*np.mean(a*(2*(X>0)-1)*(np.abs(X)>c)*(y-a*psi(X,c)-b))  ] # derivatives with respect to a,b and c
    a += -step_size*grad[0]
    b += -step_size*grad[1]
    c += -step_size*grad[2]
    c = max(0, c)
print(a,b,c)

Result : 4.823293596477329 1.9999558580638257 0.5183479531914844
[1]: https://en.wikipedia.org/wiki/Gradient_descent
