Minimize the best approximation to points that lie over a straight line with outliers for L^p norm.

line_approx.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Minimize the best approximation to points that lie over
a straight line with outliers for L^p norm.

@author: Nicolas-Guarin-Zapata
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize


def fun(a, x, y, order):
    a0, a1 = a
    return np.linalg.norm(y - a0 - a1*x, ord=order)


cm = plt.cm.get_cmap("viridis")
npts = 100
a0 = 0.0
a1 = 1.0
x = np.linspace(-1, 1, npts)
y = a0 + a1*x
y[10::20] = 10*x[10::20]
plt.plot(x, y, linewidth=0, color="black", marker=".")


orders = [0.1, 1, 2, 5, np.inf]
for cont, order in enumerate(orders):
    color = cm(cont/4)
    res = minimize(fun, [-1.0, 2.0], args=(x, y, order), method="Nelder-Mead")
    plt.plot(x, res.x[0] + res.x[1]*x, label="{}".format(order),
             color=color)
    print("{:g}: {:.4f} {:.4f}".format(order, *res.x))

plt.legend(ncol=2)
plt.show()