# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Regularied Linear Regression and Bias vs Variance

# ## 1 Import data set

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.optimize import minimize

data = loadmat('ex5data1.mat')
data


def getData(data, *args):
    return [data[i] for i in args]



X, y = data['X'], data['y']
X.shape, y.shape

Xtest, ytest = data['X'], data['ytest']
Xtest.shape, ytest.shape

Xval, yval = data['Xval'], data['yval']
Xval.shape, yval.shape


def plotData(X, y):
    plt.figure(figsize=(8, 5))
    plt.ylabel('Water flowing out of the dam (y)')
    plt.xlabel('Change in water level (X)')
    plt.plot(X, y, 'rx')
    plt.grid(True)
plotData(X, y)

X = np.insert(X, 0, 1, axis=1)
Xtest = np.insert(Xtest, 0, 1, axis=1)
Xval = np.insert(Xval, 0, 1, axis=1)


# ## 2 Regularized linear regression cost function

def h(theta, X):
    return X @ theta
def cost(theta, X, y, lam=0.):
    theta = theta.reshape((X.shape[1], 1))
    m = X.shape[0]
    a1 = h(theta, X)
    #print(a1.shape, y.shape)
    J = np.sum((a1 - y) ** 2) / (2 * m)
    #print(J)
    reg = lam * np.sum(theta[1:] * theta[1:]) / (2 * m)
    #print(reg)
    return J + reg
#theta = np.ones(2)
#cost(theta.reshape(2,1), X, y, 0.)



# ## 3 Gradient

def gradient(theta, X, y, lam=0.):
    theta = theta.reshape((theta.shape[0], 1))
    m = X.shape[0]
    grad = (X.T @ (h(theta, X) - y)) / m
    reg = float(lam) * theta / m
    reg[0] = 0.
    return (grad + reg).ravel()
#gradient(theta, X, y, 1.)


# ## 4 Training

def train(theta, X, y, lam=0.):
    fmin = minimize(cost, x0=theta, args=(X, y, lam),
                   method='TNC', jac=gradient, options={'disp':True})
    return fmin



theta1 = np.ones(2)
theta2 = np.random.rand(2)
lam = 0.2
fmin = train(theta2, X, y, lam)
print(fmin)
cost(fmin.x, X, y, lam)

plotData(X[:,1], y)
plt.plot(X[:,1], h(fmin.x, X).ravel())


# ## 5 Learning curve

def plotLearningCureve(X, y, Xval, yval, lam=0.):
    tr_cost, cv_cost = [], []
    m = list(range(1, X.shape[0] + 1))
    theta = np.random.rand(X.shape[1])
    for i in m:
        fmin = train(theta, X[:i], y[:i], lam)
        tr_cost.append(cost(fmin.x, X[:i], y[:i], lam))
        cv_cost.append(cost(fmin.x, Xval, yval, lam))
    plt.plot(m, tr_cost, label='training cost')
    plt.plot(m, cv_cost, label='cv cost')
    plt.legend(loc=1)
    plt.show()
plotLearningCureve(X, y, Xval, yval, lam=0.02)


# 模型欠拟合

# ## 6 Build poly features

def buildPolyFeatures(*args, n):
    def slove(X, n):
        m = X.shape[0]
        for i in range(2, n + 1):
            X = np.append(X, (X[:,0] ** i).reshape((m, 1)), axis=1)
        return np.insert(X, 0, np.ones(m), axis=1)
    return [slove(i, n) for i in args]
X, Xtest, Xval = buildPolyFeatures(*getData(data, 'X', 'Xtest', 'Xval'), n = 8)

X.shape



# ## 7 Standardlization

def standardlize(*args):
    def slove(X):
        for i in range(1, X.shape[1]):
            X[:,i] = (X[:,i] - X[:,i].mean()) / X[:,i].std()
        return X
    return [slove(i) for i in args]
X, Xtest, Xval = standardlize(X, Xtest, Xval)


# ## 8 Plot poly learning curve

def plotPolyLearningCurve(X, y, Xval, yval, lam=0.):
    tr_cost, cv_cost = [], []
    m = [i for i in range(1, X.shape[0] + 1)]
    theta = np.random.rand(X.shape[1])
    for i in m:
        fmin = train(theta, X[:i], y[:i], lam)
        tr_cost.append(cost(fmin.x, X[:i], y[:i]))
        cv_cost.append(cost(fmin.x, Xval, yval))
    plt.plot(m, tr_cost, label='training cost')
    plt.plot(m, cv_cost, label='cv cost')
    plt.legend(loc=1)
    plt.show()
plotPolyLearningCurve(X, y, Xval, yval)

# tr_cost too low, over fitting

plotPolyLearningCurve(X, y, Xval, yval, 1)

# ighten the fitting

plotPolyLearningCurve(X, y, Xval, yval, 100)


# less fitting

# ## 9 Find the best lam

# +
def showCost(X, y, Xval, yval, lam_vals):
    tr_cost, cv_cost = [], []
    theta = np.random.rand(X.shape[1])
    for lam in lam_vals:
        fmin = train(theta, X, y, lam)
        tr_cost.append(cost(fmin.x, X, y))
        cv_cost.append(cost(fmin.x, Xval, yval))
    print(lam_vals[np.argmin(cv_cost)])
    plt.plot(lam_vals, tr_cost, label='training')
    plt.plot(lam_vals, cv_cost, label='cross validation')
    plt.legend(loc=2)
    plt.xlabel('lambda')
    plt.ylabel('cost')
    plt.show()

lam_vals = [0, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10]
showCost(X, y, Xval, yval, lam_vals)
# -

for lam in lam_vals:
    theta = np.random.rand(X.shape[1])
    fmin = train(theta, X, y, lam)
    J = cost(fmin.x, Xtest, ytest)
    print(f'test cost lam={lam} cost={J}')


in

StarcoderdataJupyterScriptsDedupFiltered