simple liner regression with numpy

simple_liner_regression.py_with_numpy.py

import numpy as np
import matplotlib.pyplot as plt


def generate(a, b, c, n=1000):
    print("傾き {}, 切片 {}, 誤差標準偏差 {}".format(a, b, c))
    x = np.arange(n)
    y = a * x + b + c * np.random.normal(size=n)
    return x, y


def regression(x, y):
    cov = np.cov(x, y)
    Sxx, Sxy, Syy = cov[0, 0], cov[0, 1], cov[1, 1]
    # assert Sxx == np.var(x, ddof=1)  # same value, unbiased variance
    # assert Syy == np.var(y, ddof=1)  # same value, unbiased variance
    a, b, c = Sxy / Sxx, np.mean(y) - Sxy / Sxx * np.mean(x), np.sqrt(Syy - Sxy ** 2 / Sxx)
    print("傾き（推定） {}, 切片（推定） {}, 残差標準偏差 {}".format(a, b, c))
    return a, b, c


def main():
    fig = plt.figure()
    ax = fig.add_subplot(111)
    x, y = generate(0.25, -4.5, 10)
    ax.plot(x, y, label="generated")
    a, b, c = regression(x, y)
    ax.plot(x, a * x + b, label="regression")
    ax.plot(x, a * x + b + c * np.sin(x / len(x) * 100), label="regression+$1\sigma$")
    ax.legend()
    plt.show()


if __name__ == "__main__":
    main()