Linear fit when distribution of errors is highly skewed I have some datasets where the distribution of errors is expected to be highly skewed.  I'd like to do a linear fit that takes this into account.
Here is some synthetic data that shows this:
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt

n = 1000
slope = -30
x = np.random.normal(loc=0, scale=2, size=n)
y1 = slope * x + np.random.normal(loc=0, scale=100, size=n)
y2 = slope * x + scipy.stats.skewnorm.rvs(a=5, loc=0, scale=100, size=n)

plt.figure(figsize=(10,10))
plt.scatter(x,y1, alpha=0.5)
plt.scatter(x,y2, alpha=0.5)


My real data looks more similar to the orange points, not the blue points.  There's also a lot less of it, in some cases only 50-100 datapoints.  Of course I don't know that the distribution of errors is specifically skew normal - it could be better modelled by something else, eg exponentially modified Gaussian.
This is what the errors look like in the synthetic example above:

It seems when fitting a linear model to data like this, in some sense values to the left are more informative and should have a higher weight since the variance of the left half of the distribution is less.
Is there a way to do essentially a linear fit that assumes not a normal distribution of errors but some other distribution - in particular a highly skewed one?
Is it possible to simultaneously fit y=a*x+b and recover the shape of the error distribution?

There are 2 cases:

*

*If you do not know the type probability density distribution of the errors (i.e., the orange density plot), as you say, then I'd recommend to use a robust fitting procedure, such as Theil-Sen or RANSAC, rather than ordinary Least Squares (OLS). A robust estimator will reduce the chance that your estimation will go astray due to the unknown and skewed tails of the errors.


*However, if you do know the type of probability density distribution of the errors (i.e., the orange density plot), then theoretically you can do much better by using Maximum Likelihood. Assume the errors are distributed according to a known probability density $f(\epsilon; \theta)$, where $\epsilon$ are the error values and $\theta$ are a set of parameters of the family of the distribution $f$. What you need to do is mathematically express the errors and parameters as a linear regression model (as implied by the scatter-plot), and estimate all unknown parameters using Maximum Likelihood. To clarify this, see the example below.
Example for clarifying the 2nd case:
Let's assume the known type of probability density distribution of the errors is an exponential distribution, which is obviously skewed in one direction, so: $$f(\epsilon; \theta) = \theta \exp({-\theta \cdot \epsilon}) \cdot U(\epsilon)$$
where $U(\epsilon)$ is the Heavyside step function, which is zero for negative $\epsilon$ and 1 otherwise.
But we need to somehow incorporate the linear regression model implied by the scatter-plot, so for this example, I've chosen to assume: $y = a \cdot x + b + \epsilon$, where $(x, y)$ are the horizontal and vertical coordinates of a specific data-point. The error term $\epsilon$ between the $y$-axis of a data-point and its linear model is then: $\epsilon = y - (a \cdot x + b) = y - ax - b$.
Plugging this into the probability density distribution of the error for a single data-point at $(x, y)$, we have the likelihood: $$f(x, y; \theta, a, b) = \theta \exp({-\theta \cdot (y - ax - b)}) \cdot U(y - ax - b)$$ $$= \theta \exp({-\theta y + \theta ax + \theta b}) \cdot U(y - ax - b)$$
where all unknown parameters which we'll need to estimate are: $\theta, a, b$.
Let's further assume the errors $\epsilon_n$ are statistically independent for each point $(x_n, y_n)$ of the $N$ data-points, so the joint probability density of the complete dataset $(\mathbf{x}, \mathbf{y})$ of $N$ data-points is: $$f(\mathbf{x}, \mathbf{y}; \theta, a, b) = \prod_{n=1..N}{f(x_n, y_n; \theta, a, b)}$$
In order to prodceed, it would be extremely convenient to take the (natural) logarithm, so we'll have a log-likelihood function. The log-likelihood function is differentiated with respect to the unknown parameters $\theta, a, b$ in order to maximize it for the given data-points.
An iterative numerical procedure is necessary for maximization in this specific example; perhaps it will also be necessary to run it from different initializations in order to guarantee that the discovered maximum is the global one.
