线性回归练习--李宏毅机器学习hw1

2019-02-26

以前总说安全是不产生价值的岗位，结果现在又往一样不产生价值的数据分析里跳，谁让这玩意火呢。最近看机器学习，看的也不快，做个题熟悉一下，还是有不了解的东西，虽然只是个最简单的线性回归。
题目大概就是用18个指标预测第10小时pm2.5的值，介绍在https://ntumlta.github.io/2017fall-ml-hw1/
首先想的就是最简单的，只用pm2.5一个特征来预测，每天的分开算，这样就是从训练集每天取九个特征加一个答案有24-9=15种取法，一个月20天，12个月，一共122015=3600行，9列

train_data = pd.read_csv('train.csv',encoding = 'BIG5')
pm25list = train_data[train_data['測項']=='PM2.5'].iloc[:,3:]
x_tmp = []
y_tmp = []
for i in range(15):
    tempx = pm25list.iloc[:,i:i+9]
    tempx.columns=np.array(range(9))
    x_tmp.append(tempx)
    y_tmp.append(pm25list.iloc[:,i+9])
x = pd.concat(x_tmp)
y = pd.concat(y_tmp)
x_data = np.array(x,float)
y_data = np.array(y,float)

训练数据找好了，写个线性回归的函数，先用正规解算着，一会再跟梯度下降的对比。

class LinearRegression1():
    def __init__(self):
        self.coef_ = None
        self.interception_ = None
        self._theta = None
        
    def fit_normal(self,X_train,Y_train):
        X_b = np.c_[np.ones((X_train.shape[0],1)),X_train]
        self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(Y_train)
        self.coef_ = self._theta[1:]
        self.interception_ =self._theta[0]
        return self
    def predict(self,X_predict):
        assert self.coef_ is not None and self.interception_ is not None,\
        "must fit before predict"
        X_b = np.c_[np.ones((X_predict.shape[0],1)),X_predict]
        return X_b.dot(self._theta)
       
    def __repr__(self):
        return "LinearRegression1()"

代入模型求解，系数向量有9个特征加上一个b，一共十维

1
2
3

pmlrg = LinearRegression1()
pmlrg.fit_normal(x_data,y_data)
pmlrg._theta

题目有kaggle链接，虽然不能参加但是能看到得分，用的是均方根误差，定义一个

def RMSE(predict, answer):

    rmse = np.sqrt(np.sum((predict - answer)**2) /predict.shape[0])
    return rmse

计算一下和标准答案之间的RMSE
6.661983418013117
这么高，看别人都是5点多。想怎么改进一下，最简单的想法就是__加大力度__
把18个指标全都用上，跨天取值。这样读数据就比较麻烦了，因为每个月只有20天的数据，所以不能跨月。这样就有12个月，每个月一共有20天*24小时。特征依然是九个小时，取连续9小时的话最后九个小时取不了，所以每月有20*24-9=471种，12个月有12*471=5652种。
然后是列向量，只用pm2.5一个指标的时候是9个特征，现在18个指标就是18*9=162个特征。
所以x是一个5652*162的矩阵，代表一共取了5652个[九小时]，这个[九小时]里每个小时的每个污染指标都作为一个特征
y就是5652个[第十小时]

rawdata = pd.read_csv('train.csv',encoding = 'BIG5').iloc[:,3:]
data = [[] for x in range(18)]
for i in range(rawdata.shape[0]):
    if i%18 == 10:
        for d in [float(0) if x=='NR' else x for x in rawdata.iloc[i] ]:
            data[i%18].append(float(d))
    else:
        for d in rawdata.iloc[i]:
            data[i%18].append(float(d))
x = [[] for i in range(12*471)]
y = []
# 每 12 個月
for i in range(12):
    # 一個月取連續10小時的data可以有471筆
    for j in range(471):
        # 18種污染物
        for t in range(18):
            # 連續9小時
            for s in range(9):
                x[471*i+j].append(data[t][480*i+j+s] )
        y.append(data[9][480*i+j+9])
x = np.array(x)
y = np.array(y)

训练

1 2	reg = LinearRegression1() reg.fit_normal(x,y)

之后读取测试集，也不能简单的只读pm2.5的行了，一共240天，每天都是18*9=162，和训练的模型一个维度

testdata = pd.read_csv('test.csv',encoding = 'BIG5',header=None).iloc[:,2:]
xte = [[] for i in range(12*20)]
for i in range(testdata.shape[0]):
    if i%18 == 10:
        for d in [float(0) if x=='NR' else x for x in testdata.iloc[i] ]:
            xte[i//18].append(d)
    else:
        for d in testdata.iloc[i]:
            xte[i//18].append(d)

allpredict = reg.predict(xte)
rmseall = RMSE(allpredict , answer)
rmseall

6.570013915501879
尼玛就进步了这么一点点，大力也没出奇迹啊。