李宏毅机器学习HW1——PM2.5 Regression

思路主要参考：https://www.cnblogs.com/HL-space/p/10676637.html

一、作业描述

数据集描述：

1. 取12个月每个月前20天共240天的气象数据
2. 每天的数据有18个指标

二、思路分析

2.1 数据预处理

利用pd.read_csv 读取标题行之外的数据

数据中共有很多NR数据，该指标RAINFALL 为是否降雨，将NR替换为0

如图所示，取每天从0时开始到8时的所有数据为训练数据，第9个小时的PM2.5值为label，以此类推，则每天共有15个label，240天共有240X15 = 3600个label

#数据预处理，取出每天24h的数据
def dataProcessing(df): 
    df = df.iloc[:,3:].replace(['NR'],[0.0]) #将空数据替换为0    
    array = np.array(df).astype(float) # astype转换array中的元素类型
    x_list = []
    y_list = []
    for i in range(0,4320,18):
        for j in range(0,15):
            mat = array[i:i+18,j:j+9]  #每一天中每10个小时的前9个小时当作参数
            label = array[i+9,j+9] #第10列是pm2.5的值，第10个小时的pm2.5作为label
            x_list.append(mat)   #插入参数和label列表
            y_list.append(label)
    x = np.array(x_list)
    y = np.array(y_list)
    return x, y

2.2 模型建立

2.2.1 回归模型

简单起见，仅利用线性函数考虑前9个小时的PM2.5的对第10小时PM2.5的影响

其中
为第i小时的PM2.5数据，为其对应的权重值， 为偏移量， 为第10小时的预测值

2.2.2 损失函数

利用预测值和每个label的平均欧氏距离来衡量损失

其中，为label的值，为函数的预测值，为label的数量，为label的损失

未加入正则项

2.2.3 梯度下降

要想使loss最小，则需要求loss在b和w上的偏微分


更新参数

为w更新时的学习率，为b更新时的学习率

2.2.4 更新学习率


为累计梯度平方和的二次平方根，为更新前的学习率，为更新后的学习率

#调整参数，训练模型
def train(x_train,y_train,iteration):
    # ydata = b + w0 * xdata0 +...+ w8 * xdata8
    b = 0 #initial b
    w = np.ones(9) #initial w
    lr = 5 #learning rate
    pre_gra_b = 0 # adagrad history b
    pre_gra_w = np.zeros(9) # adagrad history w
    ly = len(y_train)
    
    #迭代开始, 用前9个小时的pm2.5数据计算预测第10个小时
    for i in range(iteration):
        b_grad = 0.0
        w_grad = np.zeros(9)
        
        #利用adagrad更新学习率
        for n in range(ly):
            b_grad = b_grad - 2.0 * (y_train[n] - w.dot(x_train[n,9,:])-b)
            for m in range(9):
                w_grad[m] = w_grad[m] - 2.0 * (y_train[n] - w.dot(x_train[n,9,:]) - b) * (x_train[n,9,m]) 
        
        #mean
        b_grad = b_grad/ly
        w_grad = w_grad/ly
        
        #adagrad
        pre_gra_b = pre_gra_b + b_grad ** 2
        pre_gra_w = pre_gra_w + w_grad ** 2
        
        #update parameters 
        b = b - (lr * b_grad) / np.sqrt(pre_gra_b)
        for k in w:
            k = k - (lr * w_grad ) / np.sqrt(pre_gra_w)
            
        #loss on training data, per 200
        if i % 200 == 0:
            loss = 0
            for j in range(ly):
                loss = loss + (y_train[j] - w.dot(x_train[j,9,:]) - b) ** 2
            print('iteration:{},loss:'.format(i),loss/ly)
        
    return w,b

2.3 在测试集上测试

#testing on test data         
def test(x_test,y_test,w,b):
    lt = len(y_test)
    loss = 0
    for j in range(lt):
        loss = loss + (y_test[j] - w.dot(x_test[j,9,:]) - b) ** 2
    print('the loss on test data is:', loss/lt)
    return loss/lt

三、完整代码

# -*- coding: utf-8 -*-
"""
Created on Fri Nov  1 19:49:52 2019

@author: Lenovo
"""

import pandas as pd 
import numpy as np

#数据预处理，取出每天24h的数据
def dataProcessing(df): 
    df = df.iloc[:,3:].replace(['NR'],[0.0]) #将空数据替换为0    
    array = np.array(df).astype(float) # astype转换array中的元素类型
    x_list = []
    y_list = []
    for i in range(0,4320,18):
        for j in range(0,15):
            mat = array[i:i+18,j:j+9]  #每一天中每10个小时的前9个小时当作参数
            label = array[i+9,j+9] #第10列是pm2.5的值，第10个小时的pm2.5作为label
            x_list.append(mat)   #插入参数和label列表
            y_list.append(label)
    x = np.array(x_list)
    y = np.array(y_list)
    return x, y

#调整参数，训练模型
def train(x_train,y_train,iteration):
    # ydata = b + w0 * xdata0 +...+ w8 * xdata8
    b = 0 #initial b
    w = np.ones(9) #initial w
    lr = 5 #learning rate
    pre_gra_b = 0 # adagrad history b
    pre_gra_w = np.zeros(9) # adagrad history w
    ly = len(y_train)
    
    #迭代开始, 用前9个小时的pm2.5数据计算预测第10个小时
    for i in range(iteration):
        b_grad = 0.0
        w_grad = np.zeros(9)
        
        #利用adagrad更新学习率
        for n in range(ly):
            b_grad = b_grad - 2.0 * (y_train[n] - w.dot(x_train[n,9,:])-b)
            for m in range(9):
                w_grad[m] = w_grad[m] - 2.0 * (y_train[n] - w.dot(x_train[n,9,:]) - b) * (x_train[n,9,m]) 
        
        #mean
        b_grad = b_grad/ly
        w_grad = w_grad/ly
        
        #adagrad
        pre_gra_b = pre_gra_b + b_grad ** 2
        pre_gra_w = pre_gra_w + w_grad ** 2
        
        #update parameters 
        b = b - (lr * b_grad) / np.sqrt(pre_gra_b)
        for k in w:
            k = k - (lr * w_grad ) / np.sqrt(pre_gra_w)
            
        #loss on training data, per 200
        if i % 200 == 0:
            loss = 0
            for j in range(ly):
                loss = loss + (y_train[j] - w.dot(x_train[j,9,:]) - b) ** 2
            print('iteration:{},loss:'.format(i),loss/ly)
        
    return w,b
            
#testing on test data         
def test(x_test,y_test,w,b):
    lt = len(y_test)
    loss = 0
    for j in range(lt):
        loss = loss + (y_test[j] - w.dot(x_test[j,9,:]) - b) ** 2
    print('the loss on test data is:', loss/lt)
    return loss/lt
    
def main(): 
    file = 'HW1-train.csv'
    data = pd.read_csv(file, encoding = 'gb18030')# open data
    iteration = 2000  #迭代次数
    x,y= dataProcessing(data)
    x_train, y_train = x[:3200], y[:3200]
    x_test,y_test = x[3200:3600], y[3200:3600]
    weights,bias = train(x_train, y_train,iteration)
    test(x_test,y_test,weights,bias)
    
if __name__ == '__main__':
    main()

mark：不知道为什么损失总是特别大，还需要继续调一下

四、待改进

1. 只是按天进行每十个小时的切割，应该跨天更合理一点
2. 增加正则项，防止过拟合
3. 考虑其他天气因素的影响，增加参数
4. 随机抽取训练集和测试集，而不是简单粗暴的直接前后划分