plot_lowess.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from ggplot import meat

def lowess(x, y, f=2/3., iter=3):
    n = len(x)
    r = int(np.ceil(f * n))
    h = [np.sort(abs(x - x[i]))[r] for i in range(n)]
    w = np.clip(abs(([x] - np.transpose([x])) / h), 0.0, 1.0)
    w = 1 - w * w * w
    w = w * w * w
    yest = np.zeros(n)
    delta = np.ones(n)
    for iteration in range(iter):
        for i in xrange(n):
            weights = delta * w[:, i]
            weights_mul_x = weights * x
            b1 = np.dot(weights, y)
            b2 = np.dot(weights_mul_x, y)
            A11 = np.sum(weights)
            A12 = np.sum(weights_mul_x)
            A21 = A12
            A22 = np.dot(weights_mul_x, x)
            determinant = A11 * A22 - A12 * A21
            beta1 = (A22 * b1 - A12 * b2) / determinant
            beta2 = (A11 * b2 - A21 * b1) / determinant
            yest[i] = beta1 + beta2 * x[i]
        residuals = y - yest
        s = np.median(np.abs(residuals))
        delta[:] = np.clip(residuals / (6 * s), -1, 1)
        delta[:] = 1 - delta * delta
        delta[:] = delta * delta
    return yest


# x = np.array([4,  4,  7,  7,  8,  9, 10, 10, 10, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20, 20, 22, 23, 24, 24, 24, 24, 25], np.float)
# y = np.array([2, 10,  4, 22, 16, 10, 18, 26, 34, 17, 28, 14, 20, 24, 28, 26, 34, 34, 46, 26, 36, 60, 80, 20, 26, 54, 32, 40, 32, 40, 50, 42, 56, 76, 84, 36, 46, 68, 32, 48, 52, 56, 64, 66, 54, 70, 92, 93, 120, 85], np.float)
# y_fit = lowess(x, y)

x = meat.beef.values
y = meat.veal.values
y_fit = lowess(x, y)

plt.plot(x, y, 'k.', markersize=5, label="True Data")
plt.plot(x, y_fit, '-', color='steelblue', lw=3, markersize=5, label="Lowess Fit")
plt.legend()
plt.show(1)