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

# + id="w4aUIfZjRSAe" colab_type="code" colab={}
# %matplotlib inline

# + id="LzZuMm-aTYvh" colab_type="code" colab={}
import numpy as np
import numpy.linalg as npl
import pandas as pd

import sklearn.linear_model

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')


# + [markdown] id="WktR6j1_rw3M" colab_type="text"
# # Noisy polynomials

# + id="qUwtqPO9r56o" colab_type="code" colab={}
def gen_data(n=10, sigma_noise=.2, sorted=True, fun='sin'):
    X = 2 * np.random.rand(1, n) - 1.0
    eps = np.random.randn(1, n) * sigma_noise
    if fun=='poly1':
        Y = 0.3 + .7*X + eps
    elif fun=='poly2':
        Y = 0.3 + .7*X -1.1*X**2 + eps
    elif fun =='log':
        Y = np.log(X+1)**.5 + eps
    elif fun == 'sin':
        Y = np.sin(2*np.pi*X) + eps
    else:
        raise ValueError(f"Unknown fun type: {fun}.")
    if sorted:
        xind = np.argsort(X[0, :])
        X= X[:, xind]
        Y= Y[:, xind]
    return X, Y



# + id="-CnfmEui0zu3" colab_type="code" outputId="7fafe1ab-cdaf-4146-c3aa-eda78ec5e742" executionInfo={"status": "ok", "timestamp": 1574930308906, "user_tz": -60, "elapsed": 3507, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 281}
XC,YC = gen_data(n=1000, sigma_noise=0)
fig = plt.figure()
plt.plot(XC.T,YC.T,'-') #note: plot assumes that each column is a data series
plt.title('The clean X-Y relation')
plt.show()


# + id="RrFuRqJG03ZF" colab_type="code" colab={}
def expandX(X, deg=1):
    return X**np.arange(deg + 1)[:, np.newaxis]

def expandRBF(X, N=100, sigma=1):
    Xmin = X.min()
    Xmax = X.max()
    return np.exp(-((np.linspace(Xmin, Xmax, N)[:, np.newaxis] - X) / sigma)**2)

def polyfit(X, Y, C=0.0):
    theta = npl.inv(X.dot(X.T) + np.eye(X.shape[0])*C).dot(X).dot(Y.T)
    return theta


# + id="7iib0Lst1mW7" colab_type="code" outputId="09bfb1d2-e194-4728-f652-8bb055451e19" executionInfo={"status": "ok", "timestamp": 1574930309507, "user_tz": -60, "elapsed": 4031, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 265}
plt.figure()
plt.plot(XC.T,YC.T,'-r', linewidth=2, label='clean')
for deg in [0,1,2,5,7,10,20]:
    XE = expandX(XC,deg)
    #XE = expandRBF(XC,deg)
    theta = polyfit(XE, YC)
    plt.plot(XC.T, theta.T.dot(XE).T,'-', label='deg %d' % deg, )
plt.legend(loc='upper right')
plt.show()

# + id="OyFVjPcY1mMC" colab_type="code" outputId="23c2a9af-c153-4100-9ac8-76aecee1bbe2" executionInfo={"status": "ok", "timestamp": 1574930310653, "user_tz": -60, "elapsed": 5143, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 269}
num = 20
lambda_ = 0

X,Y = gen_data(num, sigma_noise=0.4)
plt.figure()
plt.plot(X.T,Y.T,'o', label='data', alpha=1.0)
plt.plot(XC.T,YC.T,'-r', linewidth=2, label='clean')

for deg in [1, 2, 5, 10, 20, 30]:
    XE = expandX(X, deg)
    theta = polyfit(XE, Y, lambda_)
    plt.plot(XC.T, theta.T.dot(expandX(XC,deg)).T,'-', label='deg %d' % deg, )
plt.legend(loc='upper right', ncol=3)
plt.xlim([-1.1, 1.1])
plt.ylim([-2,2.5])
plt.show()

# + id="FRRVrVJu1lb_" colab_type="code" outputId="bea879a4-ff8a-4339-acf8-4c7f34cabd4e" executionInfo={"status": "ok", "timestamp": 1574930310654, "user_tz": -60, "elapsed": 5055, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 281}
#Hint: change the degree and number of samples
deg = 50
num = 15
lambda_ = 1e-2

plt.figure()
plt.title('Various degree=%d fits' % deg)
plt.plot(XC.T, YC.T, '-r', linewidth=2, label='Clean rel.')
for run in range(5):
    X,Y = gen_data(num, sigma_noise=0.4)
    XE = expandX(X,deg)
    theta = polyfit(XE, Y, lambda_)
    XCE = expandX(XC,deg)
    plt.plot(XC.T, theta.T.dot(XCE).T,'-' )

plt.xlim([-1.1, 1.1])
plt.ylim([-2, 2])
plt.show()

# + [markdown] id="GszVgLYHr1zY" colab_type="text"
# # Outliers and robust regression

# + id="O2DwA8ka6Ai5" colab_type="code" outputId="e498c8de-099f-4459-8947-1e0aec6432a9" executionInfo={"status": "ok", "timestamp": 1574930310664, "user_tz": -60, "elapsed": 5026, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 34}
XC.shape

# + id="kKblHNkO5FGh" colab_type="code" outputId="a52e94c1-32c4-4279-ab7b-68a21f3b8dd3" executionInfo={"status": "ok", "timestamp": 1574930310669, "user_tz": -60, "elapsed": 4999, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 281}
XC,YC = gen_data(n=30, sigma_noise=0.2, fun='poly1')
# now append some outliers

XO = np.concatenate([XC, [[-1, -0.9, -0.75]]], axis=1)
YO = np.concatenate([YC, [[3, 3, 3]]], axis=1)
fig = plt.figure()
plt.plot(XO.T,YO.T,'o') #note: plot assumes that each column is a data series
plt.title('X-Y relation with outliers')
plt.show()

# + id="1VUFhrjS6fpg" colab_type="code" outputId="84b79ab6-03bc-4d63-cc48-5705abd88b50" executionInfo={"status": "ok", "timestamp": 1574930311163, "user_tz": -60, "elapsed": 5070, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 353}
fig = plt.figure()
plt.plot(XO.T,YO.T,'o') #note: plot assumes that each column is a data series

XE = expandX(XC, 1)
poly = sklearn.linear_model.LinearRegression()
YFit = poly.fit(XE.T, YC.T)
plt.plot(XC.T, poly.predict(XE.T), label='Lin Reg no outliers')

poly = sklearn.linear_model.HuberRegressor()
YFit = poly.fit(XE.T, YC.T)
plt.plot(XC.T, poly.predict(XE.T), label='Huber Reg no outliers')

# Models on outliers
if 1:
    XE = expandX(XO, 1)
    poly = sklearn.linear_model.LinearRegression()
    YFit = poly.fit(XE.T, YO.T)
    plt.plot(XO.T, poly.predict(XE.T), label='Lin Reg with outliers')

    poly = sklearn.linear_model.HuberRegressor()
    YFit = poly.fit(XE.T, YO.T)
    plt.plot(XO.T, poly.predict(XE.T), label='Huber Reg with outliers')

_ = plt.legend()

# + [markdown] id="YjVuX0bxrzuQ" colab_type="text"
# # Forward stagewise

# + id="T5utaWkdTjn-" colab_type="code" outputId="a1bd6dbe-86e5-4c79-9593-dbac4f0f3d31" executionInfo={"status": "ok", "timestamp": 1574930312245, "user_tz": -60, "elapsed": 6110, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 238}
# Load the data

cancer_df = pd.read_csv(
    'https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data',
    sep='\t', header=0, index_col=0)

train_id = cancer_df['train']
del cancer_df['train']

# Standardize the whole data, as per its readme
# https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.info.txt
# NB, usually you should standardize based on train set only!!

cancer_df -= cancer_df.mean(0, numeric_only=True)
cancer_df /= cancer_df.std(0, numeric_only=True)

cancer_train_df = cancer_df[train_id=='T']
cancer_test_df = cancer_df[train_id=='F']

X_cols = [c for c in cancer_df.columns if c!='lpsa']
Y_col = 'lpsa'

print(f'Train nas {len(cancer_train_df)} saples\n'
      f'Test nas {len(cancer_test_df)} saples')

cancer_train_df.head()

# + id="oJSNLK8etmVN" colab_type="code" outputId="46018b5e-2294-4839-e638-2db1044074d9" executionInfo={"status": "ok", "timestamp": 1574930312247, "user_tz": -60, "elapsed": 6078, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 265}
# Forward stagewise

epsilon = 1e-1
theta = np.zeros((len(X_cols),))
X = np.array(cancer_train_df[X_cols])
r = np.array(cancer_train_df[Y_col])

coefs = [np.array(theta)]
errs = []
last_err = np.sum(r**2)

for _ in range(1000000):
    # compute correlations
    residual_correlations = r @ X
    best_feat = np.argmax(np.abs(residual_correlations))
    max_corr = residual_correlations[best_feat]
    
    theta[best_feat] += epsilon * np.sign(max_corr)
    r -= epsilon * np.sign(max_corr) * X[:, best_feat]
    cur_err = np.sum(r**2)

    if last_err - cur_err < 0:
        break
    last_err = cur_err
    errs.append(cur_err)
    coefs.append(np.array(theta))


coefs = np.array(coefs).T
tau = np.abs(coefs).sum(0)

for i, name in enumerate(X_cols):
    plt.plot(tau, coefs[i], label=name)
#_ = plt.legend()

# + id="ZmbKn55tUTUs" colab_type="code" cellView="both" outputId="1f4f669c-266e-4741-e2a9-767697e82781" executionInfo={"status": "ok", "timestamp": 1574930312248, "user_tz": -60, "elapsed": 6045, "user": {"displayName": "<NAME>", "photoUrl": "", "userId": "18100051538790809279"}} colab={"base_uri": "https://localhost:8080/", "height": 265}
#@title  Lasso LARS

(alphas, active, coefs
 ) = sklearn.linear_model.lars_path(np.array(cancer_train_df[X_cols]), 
                                    np.array(cancer_train_df[Y_col]))
 
# coefs is D x n_steps
tau = np.abs(coefs).sum(0)

for i, name in enumerate(X_cols):
    plt.plot(tau, coefs[i], label=name)
_ = plt.legend()

# + id="3YHQ81B9ZotY" colab_type="code" colab={}


StarcoderdataJupyterScriptsDedupFiltered