<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Cost function for heavy-tailed noise
Step1:  Generate some data
Our "ground truth" data is a simple linear function.
Step2:  Let's try adding Gaussian (normal) noise...
Step3:  ... or some Cauchy (heavy-tailed) noise
Step4:  Trying to recover the linear function using Ridge regression
Gaussian noise
For Gaussian noise the recovery works fairly well.
Step5:  Cauchy noise
Let's try the same with the Cauchy noise.
Step6:  That looks much less impressive. The problem is that the squared-$\ell^2$-norm minimized by Ridge regression is too sensitive to outliers.
Using a less sensitive cost function
Let
$$z_i = (X \cdot W - y)_i$$
denote the residuals. Instead of the usual least-squares approach of minimizing
$$\min_W ||z||_{\ell^2}^2 = \frac12 \sum_i |z_i|^2 \,, $$
let's try to minimize the $\ell^1$ norm
$$\min_W ||z||_{\ell^1} = \sum_i |z_i|$$
instead. Alternatively, since that is not smooth, we could try to minimize an objective like
$$ \min_W \eta^2 \sum_i \left( \sqrt{ 1 + \left( \frac{|z_i|}{\eta} \right)^2 } - 1 \right) \,, $$
where $\eta$ is a positive scale parameter.
This is a smooth convex cost function based component-wise on the Japanese bracket
Step7:  $\ell^1$ cost
If we switch to the $\ell^1$ cost, we start to get very good results!
Step8:  Japanese cost
If we switch to the Japanese bracket cost, we get results just as good. For this simple example, the $\ell^1$ cost already gave us good answers, but for more complex problems the smoothness of the Japanese bracket may make it a numerically safer choice.
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function
import numpy as np
from sklearn.linear_model import Ridge
from flexible_linear import FlexibleLinearRegression
import matplotlib
import matplotlib.pyplot as plt
matplotlib.style.use('ggplot')
%matplotlib inline
np.random.seed(1)
Explanation: Cost function for heavy-tailed noise
End of explanation
N = 500
A = 50
B = 3
x = np.linspace(0, 100, N, dtype=float)
X = x.reshape(-1, 1)  # scikit-learn wants 2d arrays
y = A + B * x
plt.plot(x, y, '-')
Explanation: Generate some data
Our "ground truth" data is a simple linear function.
End of explanation
y_gauss = y + 20*np.random.randn(N)
plt.plot(x, y, '-', x, y_gauss, '-')
plt.legend(['True', 'Noisy'], loc='upper left')
Explanation: Let's try adding Gaussian (normal) noise...
End of explanation
y_cauchy = y + 20*np.random.standard_cauchy(N)
plt.plot(x, y, '-', x, y_cauchy, '.')
plt.legend(['True', 'Noisy'], loc='lower right')
Explanation: ... or some Cauchy (heavy-tailed) noise
End of explanation
clf = Ridge()
clf.fit(X, y_gauss)
pred = clf.predict(X)
plt.plot(x, y, '-', x, pred, '-')
plt.legend(['True', 'Recovered'], loc='upper left')
print("     True: %.3f + %.3f * x" % (A, B))
print("Recovered: %.3f + %.3f * x" % (clf.intercept_, clf.coef_[0]))
Explanation: Trying to recover the linear function using Ridge regression
Gaussian noise
For Gaussian noise the recovery works fairly well.
End of explanation
clf = Ridge()
clf.fit(X, y_cauchy)
pred = clf.predict(X)
plt.plot(x, y, '-', x, pred, '-')
plt.legend(['True', 'Recovered'], loc='upper left')
print("     True: %.3f + %.3f * x" % (A, B))
print("Recovered: %.3f + %.3f * x" % (clf.intercept_, clf.coef_[0]))
Explanation: Cauchy noise
Let's try the same with the Cauchy noise.
End of explanation
clf = FlexibleLinearRegression(cost_func='l2', C=0.0)
clf.fit(X, y_cauchy)
pred = clf.predict(X)
plt.plot(x, y, '-', x, pred, '-')
plt.legend(['True', 'Recovered'], loc='upper left')
print("     True: %.3f + %.3f * x" % (A, B))
print("Recovered: %.3f + %.3f * x" % (clf.coef_[0], clf.coef_[1]))
Explanation: That looks much less impressive. The problem is that the squared-$\ell^2$-norm minimized by Ridge regression is too sensitive to outliers.
Using a less sensitive cost function
Let
$$z_i = (X \cdot W - y)_i$$
denote the residuals. Instead of the usual least-squares approach of minimizing
$$\min_W ||z||_{\ell^2}^2 = \frac12 \sum_i |z_i|^2 \,, $$
let's try to minimize the $\ell^1$ norm
$$\min_W ||z||_{\ell^1} = \sum_i |z_i|$$
instead. Alternatively, since that is not smooth, we could try to minimize an objective like
$$ \min_W \eta^2 \sum_i \left( \sqrt{ 1 + \left( \frac{|z_i|}{\eta} \right)^2 } - 1 \right) \,, $$
where $\eta$ is a positive scale parameter.
This is a smooth convex cost function based component-wise on the Japanese bracket:
$$ \langle z \rangle := \sqrt{1 + |z|^2} \,. $$
It behaves like the least-squares cost function for $|z| \ll \eta$,
but grows like $O(|z|)$ for $|z| \gg \eta$.
So it gives us the best of $\ell^1$ and squared-$\ell^2$:
The smoothness of squared-$\ell^2$ for small values
and the insensitivity to outliers of $\ell^1$ for large values.
We use a custom FlexibleLinearRegression scikit-learn estimator that allows us to try these different cost functions.
Note: For simplicity, we omit regularization in the above (and below by setting C=0.0). It really does not matter for this toy problem.
squared $\ell^2$ cost
This is pretty much the same as Ridge above (minus the regularization). It performs horribly.
End of explanation
clf = FlexibleLinearRegression(cost_func='l1', C=0.0)
clf.fit(X, y_cauchy)
pred = clf.predict(X)
plt.plot(x, y, '-', x, pred, '-')
plt.legend(['True', 'Recovered'], loc='upper left')
print("     True: %.3f + %.3f * x" % (A, B))
print("Recovered: %.3f + %.3f * x" % (clf.coef_[0], clf.coef_[1]))
Explanation: $\ell^1$ cost
If we switch to the $\ell^1$ cost, we start to get very good results!
End of explanation
clf = FlexibleLinearRegression(cost_func='japanese', C=0.0, cost_opts={'eta': 10.0})
clf.fit(X, y_cauchy)
pred = clf.predict(X)
plt.plot(x, y, '-', x, pred, '-')
plt.legend(['True', 'Recovered'], loc='upper left')
print("     True: %.3f + %.3f * x" % (A, B))
print("Recovered: %.3f + %.3f * x" % (clf.coef_[0], clf.coef_[1]))
Explanation: Japanese cost
If we switch to the Japanese bracket cost, we get results just as good. For this simple example, the $\ell^1$ cost already gave us good answers, but for more complex problems the smoothness of the Japanese bracket may make it a numerically safer choice.
End of explanation
<END_TASK>