%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.preprocessing import PolynomialFeatures
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import make_pipeline

sns.set()
plt.rc('font', family='SimHei')
plt.rc('axes', unicode_minus=False)

x_fit = (np.random.rand(100) * 2 * np.pi)[:, np.newaxis]
y_fit = np.sin(x_fit) + 0.2 * np.random.rand(100)[:, np.newaxis]

plt.figure(figsize=(10, 10))
plt.plot(x_fit, y_fit, 'o')

poly = PolynomialFeatures(degree=7)
model = LinearRegression(fit_intercept=True)
model.fit(poly.fit_transform(x_fit), y_fit)

print(f'斜率拟合结果：\n{model.coef_}\n')
print(f'截距拟合结果：\n{model.intercept_}\n')

res = model.predict(poly.fit_transform(np.linspace(-1, 7, 100)[:, np.newaxis]))

plt.plot(np.linspace(-1, 7, 100), res)
plt.title('使用多项式基函数拟合正弦波（线性回归）')

class GaussianFeatures(BaseEstimator, TransformerMixin):
    def __init__(self, N:int, width_factor:float = 2.0):
        self.N = N
        self.width_factor=width_factor
        self.centers_ = None
        self.width_ = None

    @staticmethod
    def _gauss_basis(x, y, width, axis=None):
        arg = (x - y) / width
        return np.exp(-0.5 * np.sum(arg ** 2, axis))

    def fit(self, x, y=None):
        self.centers_ = np.linspace(x.min(), x.max(), self.N)
        self.width_ = self.width_factor * (self.centers_[1] - self.centers_[0])
        return self

    def transform(self, x):
        return self._gauss_basis(x[:, :, np.newaxis], self.centers_, self.width_, axis=1)


model = make_pipeline(GaussianFeatures(20), LinearRegression())
x_fit = np.random.rand(100) * 2 * np.pi
y_fit = np.sin(x_fit) + 0.2 * np.random.rand(100)
x_fit = x_fit[:, np.newaxis]; y_fit = y_fit[:, np.newaxis]
x_test = np.linspace(0, 2 * np.pi, 200)[:, np.newaxis]

model.fit(x_fit, y_fit)
res = model.predict(x_test)

plt.figure(figsize=(10, 10))
plt.plot(x_fit, y_fit, 'o')
plt.plot(x_test, res)
plt.title('使用高斯基函数拟合正弦波（线性回归）')

x_fit = np.random.rand(60) * 2 * np.pi - np.pi
y_fit = np.sin(x_fit) + 0.3 * np.random.rand(60)
x_test = np.linspace(-np.pi, np.pi, 200)
x_fit, y_fit, x_test = x_fit[:, np.newaxis], y_fit[:, np.newaxis], x_test[:, np.newaxis]

model1 = make_pipeline(GaussianFeatures(30), LinearRegression())
model2 = make_pipeline(GaussianFeatures(30), Ridge(alpha=0.2))

model1.fit(x_fit, y_fit); model2.fit(x_fit, y_fit)
res1, res2 = model1.predict(x_test), model2.predict(x_test)

fig, axs = plt.subplots(4, 1, figsize=(15, 20))
ax_linear = axs[0]  # type: plt.Axes
ax_coef_linear = axs[1]  # type: plt.Axes
ax_ridge = axs[2]  # type: plt.Axes
ax_coef_ridge = axs[3]  # type: plt.Axes

ax_linear.plot(x_fit, y_fit, 'o')
ax_linear.plot(x_test, res1)
ax_linear.set_xlim(-3.5, 3.5)
ax_linear.set_ylim(-1.5, 1.5)
ax_linear.set_title('简单线性回归结果')

ax_coef_linear.plot(model1.steps[0][1].centers_[:, np.newaxis], model1.steps[1][1].coef_.reshape(30, 1))
ax_coef_linear.set_title('简单线性回归振幅')

ax_ridge.plot(x_fit, y_fit, 'o')
ax_ridge.plot(x_test, res2)
ax_ridge.set_xlim(-3.5, 3.5)
ax_ridge.set_ylim(-1.5, 1.5)
ax_ridge.set_title('岭回归结果')

ax_coef_ridge.plot(model2.steps[0][1].centers_[:, np.newaxis], model2.steps[1][1].coef_.reshape(30, 1))
ax_coef_ridge.set_title('岭回归振幅')

x_fit = np.random.rand(60) * 2 * np.pi - np.pi
y_fit = np.sin(x_fit) + 0.3 * np.random.rand(60)
x_test = np.linspace(-np.pi, np.pi, 200)
x_fit, y_fit, x_test = x_fit[:, np.newaxis], y_fit[:, np.newaxis], x_test[:, np.newaxis]

model1 = make_pipeline(GaussianFeatures(30), Ridge(0.2))
model2 = make_pipeline(GaussianFeatures(30), Lasso(0.002))

model1.fit(x_fit, y_fit); model2.fit(x_fit, y_fit)
res1, res2 = model1.predict(x_test), model2.predict(x_test)

fig, axs = plt.subplots(4, 1, figsize=(15, 20))
ax_linear = axs[0]  # type: plt.Axes
ax_coef_linear = axs[1]  # type: plt.Axes
ax_ridge = axs[2]  # type: plt.Axes
ax_coef_ridge = axs[3]  # type: plt.Axes

ax_linear.plot(x_fit, y_fit, 'o')
ax_linear.plot(x_test, res1)
ax_linear.set_xlim(-3.5, 3.5)
ax_linear.set_ylim(-1.5, 1.5)
ax_linear.set_title('岭回归结果')

ax_coef_linear.plot(model1.steps[0][1].centers_[:, np.newaxis], model1.steps[1][1].coef_.reshape(30, 1))
ax_coef_linear.set_title('岭回归振幅')

ax_ridge.plot(x_fit, y_fit, 'o')
ax_ridge.plot(x_test, res2)
ax_ridge.set_xlim(-3.5, 3.5)
ax_ridge.set_ylim(-1.5, 1.5)
ax_ridge.set_title('Lasso回归结果')

ax_coef_ridge.plot(model2.steps[0][1].centers_[:, np.newaxis], model2.steps[1][1].coef_.reshape(30, 1))
ax_coef_ridge.set_title('Lasso回归振幅')