815惩罚回归之Boston_housing exercise

在tutorial1里，part b是关于boston housing的练习，这个练习主要有以下几个步骤：
1.读取数据，构建自变量的平方项和三次方项，为自变量ZN构建一个哑变量
2.对数据进行标准化，跑一个最小二乘回归
3.使用Lasso做惩罚回顾
4.使用交叉验证做Lasso

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

@author: Clementine

"""

# In[1]:
#首先是我们熟悉的导入包过程
#可以看到有些包用都没用，不知道他为什么要导入，但是他要导入就导入吧
import numpy as np
import pandas as pd
import seaborn as sns
import requests
from matplotlib import pyplot as plt
from urllib.request import urlopen
sns.set()
from sklearn.datasets import load_boston #使用sklearn连接boston数据库
from sklearn import linear_model #sklearn是一个机器学习包，这里导入的是线性模型

x,y=load_boston(return_X_y=True) #从数据库导入数据

y=pd.DataFrame(y, columns=['medv']) #把数据转换成pandas的dataframe类型
x=pd.DataFrame(x, columns=["crim","zn","indus","chas","nox","rm","age","dis","rad","tax","ptratio","b","lstat"])

sns.set(rc={'figure.figsize':(11.7,8.27)}) #设置图片参数
sns.distplot(y, bins=30) #生成y的图像
plt.show() #输出y的图像

没啥好说的就是把xy读了进来，然后把y的分布画了出来。

# In[2]:

x2=x**2 #x平方项
x2.columns=x2.columns+'^2'#一开始x2的列名和x的列名一样，这里把列名改成xxx^2,即crim变成crim^2

x3=x**3 #x3次方项
x3.columns=x3.columns+'^3'

dummy=(x['zn']>0)*1 #x中zn列对应的哑变量，如果zn为正，则dummy为1，否则为0
dummy.columns=['zn Dummy']

newX=pd.concat((x, dummy, x2, x3), axis=1) #把x,dummy,x2,x3合成新的dataframe （newX)
#前面改列名就是为了这里，如果一个dataframe里面有同样名称的列名那么调用的时候会出问题的。（虽然我们也没有调用，但是也要严谨

mean=np.mean(newX, axis=0) #计算每列的平均值
std=np.std(newX, axis=0) #计算每列的方差

stdX=(newX-mean)/std #对newX进行标准化

np.std(stdX.iloc[:,1])
#计算第一列的方差，是1哟~（因为标准化了，所以是均值为0方差为1

这里输出了1

# In[3]:

XX=pd.concat((pd.DataFrame(np.ones([newX.shape[0],1]), columns=['intercept']), stdX), axis=1)
#新加一列自变量，全为1，那么这一列对应的beta就是截距~哒哒~

beta=np.linalg.inv(XX.T@XX)@(XX.T@y)
#为了让我们的XX有点存在感，俺把它替换上去了，但是它对应的beta好像是0
#beta=np.linalg.inv(stdX.T@stdX)@(stdX.T@y)
plt.plot(beta)#画出beta

这里要专门解释一下，老师写的代码里面没有使用XX，都是用的stdX，俺觉得不对，俺就要用XX，当然，你也可以用stdX，这样就能得到和老师一样的答案了。
如图是OLS的beta，OLS就是使用np.linalg.inv做的，就是矩阵算出来的。

下面是使用Lasso做惩罚回归

# In[4]:

model=linear_model.Lasso(alpha=0.5, tol=1e-4, max_iter=10000)
model.fit(XX,y)

print(model.coef_)

得到如下系数
[ 0. -0. 0. -0. 0.0039127 -0.
0. -0. -0.04050466 -0. -0. -1.41279353
0.54466037 -3.32007833 0. -0.21952744 0. -0.
0.20497481 -0. 0. -0. -0. -0.
-0. -0. 0. -0. -0. 0.
-0. 0.15540629 -0.17625958 3.60267302 -0. -0.
-0. -0. -0. 0. 0. ]

然后我们使用如下语句查看哪些自变量被删掉了，即beta=0

# In[5]:

np.where(model.coef_==0)

得到的结果是
(array([ 0, 1, 2, 3, 5, 6, 7, 9, 10, 14, 16, 17, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30, 34, 35, 36, 37, 38, 39, 40],
dtype=int64),)

随后，将beta系数画出来

# In[6]:

idx=np.where(model.coef_!=0)
XX.iloc[:,list(idx[0][:])]
plt.plot(model.coef_)

然后是交叉验证的Lasso，关于交叉验证是验证什么，我个人认为是lambda的选择，我看那个函数的参数没看懂，有想要弄懂的可以到
偷来的sklearn.linear_model.LassoCV()网页
如果你看懂了，麻烦也给我讲讲，谢谢（狗头）

这里我们做的是10-fold cross validation，关于这个方法，我想我们是学过的，不过考虑到说不定有网友来看我的笔记，我还是
随便偷了个cross validation链接
方便大家理解什么是交叉验证

# In[7]:

model=linear_model.LassoCV(cv=10, random_state=5, tol=1e-10, max_iter=10000000).fit(XX,y.values.ravel())
print(model.coef_)

得到的结果是：
[ 0. -2.86520855 -2.23432575 -1.09056479 0.0469 -0.74657588
-5.4953677 -0.34864818 -4.91586789 3.08878787 -1.9710501 -3.51545216
1.52223698 -9.53535708 0.99927966 1.26852115 -0. 1.21409531
0.05863961 -1.56442394 0. -0. 2.32359789 0.
-0. 0. 0. 5.37498803 -0. 1.80580263
0. 0.55458561 -0. 7.91782261 0.40985549 0.
0. -0. 1.89813824 -1.05774706 -0. ]

定位被退化为0的beta

# In[8]:

np.where(model.coef_==0)

得到：
(array([ 0, 16, 20, 21, 23, 24, 25, 26, 28, 30, 32, 35, 36, 37, 40],
dtype=int64),)

最后一步啦~画出beta

# In[9]:

idx=np.where(model.coef_!=0)
XX.iloc[:,list(idx[0][:])]
plt.plot(model.coef_)

请自行思考这些beta图的不同说明了什么，因为我也不知道 （狗头）

我宣布，下班！