815线性惩罚回归笔记

这学期上了一门课，鉴于老师写完代码不检查以及很多同学不懂编程，我决定详细地分享每一个代码，今天就从tutorial1的part a开始。
鉴于很多同学是第一次接触python，我们把装包在最开始讲一遍
如果你跑代码出现错误：
ModuleNotFoundError: No module named ‘xxx’
你可以把它复制然后百度一下（狗头）
好吧，其实就是需要安装包，如果安装了Anaconda的同学就比较快乐了
给你们偷来的anaconda装包链接
我上面的xxx就是链接里的图片中的包名，具体就是pandas，matplotlib，tensorflow这些
如果没有安装Anaconda的同学也不用担心，也不一定每个人用cmd安都会出问题的
我也给你们偷了个用cmd装包的链接
好的，装包问题差不多就到这里结束了，如果装包后仍有什么问题，你可以尝试在各种地方找到我的联系方式，网络一线牵（狗头）

在放代码之前，简单给朋友们讲一下，老师给的代码是一段一段跑的，如果你都写在一个文件里，就需要选中一段（按住鼠标左键选，理解一下），然后选择 运行->运行选中代码或当前行
但是我重写了，就是说我的分段是按有输出就分个段，方便我把输出放上来这样子

# -*- coding: utf-8 -*-
"""
Created on Fri Jan 28 13:10:09 2022

@author: Clementine
"""

# In[1]:
#这下面一堆乱七八糟的东西，不用知道是什么，反正就是。。。包
#有兴趣可以网上冲浪了解一下到底是干什么的
import numpy as np
import scipy as sp
import pandas as pd
import pdb
import seaborn as sns
from sklearn.preprocessing import StandardScaler
#上面这个是数据预处理的包中的标准化的函数
from matplotlib import pyplot as plt #这是画图的包
sns.set()
#如果你用的pyCharm的话还需要加上：
#get_ipython().run_line_magic('config',"InlineBackend.figure_format = 'retina' # this is used to make the plots high-res, if run in PyCharm comment this out, it will give an error.")
#如果你用的IDLE或者spyder就不需要了，其他的我没用过
plt.rcParams["axes.grid"]=False
#这个大概是图的参数设定，想知道的可以搜索matplotlib.pyplot.rcParams

#首先是获得一些beta，我们设置一个函数
#注意我们设置的这个函数在使用的时候p要大于1
def get_sparse_beta(p):
    beta=np.zeros([p,1])
    #先声明一个p*1的矩阵，里面装的都是0 （zeros)
    for j in range(1,p):
        beta[j-1]=1/j**2 #从beta[0]就是beta里的第一个空间，存进去1/j^2
        #即beta0=1， beta1=1/4以此类推，所以得到了p个beta,beta p-1就很小很小。
    return beta


stand=1#这个代码开始讲废话了，反正就是要你必须归一化，不如直接搞，非要整个if判断
#但是其实呢这里是告诉你在更多数据的情况下，数据可以有一个flag表示需不需要做归一化，stand可以成为函数的输入。（不重要）
T=100 #时间，也可以看做样本数量
p=200 #beta，即自变量个数
beta=get_sparse_beta(p) #生成p个beta 注意这里的get_sparse_beta()要与之前def的名字一样
#而且def里面的p可以用别的字母替代
#就是说其实上面那个def可以写成def a(p)，那么这里写成m=200, beta=a(m)也可
#上面这段话看不懂也不重要，以后会懂的（狗头）

#下面是把beta画出来
plt.plot(beta) #先画一条线
plt.ylabel('beta') #加上y轴的名字
plt.xlabel('p') #加上x轴的名字

然后我们就得到的200个beta
（才发现CSDN会自动给图片加水印，好厉害）

我们已经有beta系数了，下一步就是生成自变量和应变量。
好，继续！冲！

# In[2]:

X=np.random.randn(T,p)#生成（大概是正态分布的）随机数，就是p个自变量，每个自变量是长达T的时间序列
e=np.random.randn(T,1)#生成1个长达T的时间序列，随机扰动项
y=X@beta+e #y=betaX+error,@大概是说要按照时间序列来，不要搞成矩阵相乘了。

if stand:#这里是一个标准化的过程，使xy正态分布
    sc=StandardScaler(copy=True, with_mean=True, with_std=False)
    y=sc.fit_transform(y)
    X=sc.fit_transform(X)

#下面是一通操作让我们发现ols在p>T的时候非常智障
bols=np.linalg.inv(X.T@X)@(X.T@y)#这个就是直接用ols计算的beta
plt.plot(beta,label='real')#画出真实的beta
plt.plot(bols)#画出计算出的beta
plt.legend(['real','ols'])#加上标签

然后我们就得到了

我们发现ols出来的系数可达上千（然而我们实际的beta系数最大的才是1），所以在p>>T的时候是不好直接用ols的哟~

在我们已经知道OLS不太聪明的时候，我们就想知道，惩罚回归到底做得怎样。
首先我们使用0范式来试试！

# In[3]:

lam=100 #可以改变lambda的值看看结果如何
bridge=np.linalg.inv(X.T@X+lam*np.eye(p))@(X.T@y)
#这里np.eye(p)是一个p*p的单位矩阵
#np.linalg.inv是求矩阵的逆
#所以这里是用0范式做的惩罚回归
plt.plot(bridge)
plt.plot(beta)

哇，much better，但是我们发现基本上所有的系数都是接近0，但是还是不好，跟实际相差还是很大

然后试试2范式

# In[4]:
#下面开始使用Ridge和Lasso做惩罚回归
from sklearn.linear_model import Ridge, Lasso
clf = Ridge(alpha=lam)
clf.fit(X,y);
plt.plot(clf.coef_[0],label='package')#这是用2范式（即Ridge）做的
plt.plot(bridge,label='closed form')#这是0范式
plt.legend(['package','closed form'])
#0范式很好地遮住了2范式的蓝色线条而已，它是在的

这次画的图是0范式和2范式，发现结果很像嘛

下面，终于到了万众瞩目的1范式，Lasso!

# In[5]:
lam=1e-2#我没办法解释为什么要在这里改变lambda，老师想变吧
lasso=Lasso(alpha=lam)
lasso.fit(X,y)
lasso.coef_
plt.plot(lasso.coef_);
plt.plot(beta)
plt.legend(['lasso','real'])

#下班！

如果把较小的beta看作0的话，可以说lasso把不少beta置零了
这次tutorial告诉我们，lasso可以起到variable selection的作用

下班！
part b下次再说