李巨集毅機器學習:PM2.5預測MATLAB實現

上一篇採用兩個引數進行gradient descend計算線性迴歸問題，由於自己也是新手才學不久，最近又在啃周志華著的《機器學習》，深深的陷進去了，現在才想起吧這個未完成的任務完成，我也是新手也請大神指點。
直接進入正題

將資料提取出來我這裡分成了兩步

1.將資料從train.csv中提取到MATLAB中

將資料存放在data矩陣中，其中矩陣的每一行是一種氣象資料（feature）。tip:MATLAB中csvreader函式，計算行和列是從0開始計數的
BTW:由於MATLAB中csvreader函式只能讀取資料，我先將train.csv中的NR改成了0（我也很痛苦，本來想直接實現的，試了好多辦法比如textscan都不太好用，到時候學習中有好辦法再改吧）

直接上程式碼

new_feature=1;%定義CSV檔案中特徵行的變化 取1是指第二行才有資料

for day=0:239%一次迴圈讀取1天24小時1種特徵  定義收集資料天數的變化，共240天
    for feature=1:18

        data( feature,day*24+1:day*24+24 )=...%讀取的是1*24列的矩陣
             csvread('train.csv',new_feature,3,[new_feature,3, new_feature,26 ]);
        new_feature=new_feature+1
 
;

    end

end
clear day feature new_feature

2.將data矩陣中的資料對映到x,y中

選取前9小時氣象資料作為預測的輸入（x），預測當前小時的PM2.5
x是5652X163階矩陣，5652是指train.csv中採集的是一個月中前20天的資料，每天24小時，取前9小時作為訓練資料共有24X20-9組資料又有12個月
同時為了方便矩陣計算將x的最後一列加上一個全為1的向量用於計算bias。
程式碼如下

%data中前24*20為一個月，對每個月進行取值每10小時1組前9為training set共471(480-9)組
x=zeros(471
 
*12,18*9);
y=zeros(471*12,1);
for month=0:11
    for j=1:471
        for hours=0:9
            if(hours<9)
                x(471*month+j,hours*18+1:hours*18+18)=data(:,480*month+j +hours);%前9組做輸入
                %480*month+j是指每個月有480組資料，每個月內的資料由j選擇
            else
                y(471*month+j,1)=data(10,480*month+j +hours);%PM2.5在第十行
            end
        end
    end
end

%將全為1的向量加入矩陣末尾用於計算bias
bias=ones(471*12,1);
x=[x,bias];

clear month j hours

固定學習速率計算

先進行矩陣計算
符號規定X為輸入矩陣（包含最後一列為1），Y∗train.csv的實際值，Ycalc為計算值，w為係數向量，最後一個向量為bias;L為損失函式則

程式碼如下

%initial parameters
w=zeros(18*9+1,1);%共19*9個weight和1個bias
%y_calcu=x*w 
lr=10e-11;
iteration=10000;

%figure(1);%用於顯示cost函式的變化
%hold on

for i=0:iteration
    detal_y =x*w-y;
    grad=x'*detal_y;
    cost= (sum(detal_y.^2/length(y)))^0.5;
    w=w-lr*grad;
%   scatter(i,cost,'r');
    fprintf('iteration=%8d\t|\tcost=%f\n',i,cost);
end

%hold off

Tip:這裡的學習速率同我上篇文章一樣是測試出來的，再大一個數量級loss function馬上就發散。
相關理論推薦看這篇文章：梯度下降（Gradient Descent）小結

Adagrad優化學習速率

%initial parameters
w=zeros(18*9+1,1);%共19*9個weight和1個bias
%y_calcu=x*w 
lr=1000;
iteration=1e4;
sum_grad=0;%保證分母不為0

% figure(1);%用於顯示cost函式的變化
% hold on

for i=0:iteration
    detal_y =x*w-y;
    grad=x'*detal_y;
    sum_grad=sum_grad+grad.^2;
    cost= (sum(detal_y.^2/length(y)))^0.5;
    w=w-lr*grad./(sum_grad).^0.5;
%     scatter(i,cost,'r');
    fprintf('iteration=%8d\t|\tcost=%f\n',i,cost);

    if(max( abs(grad) ) <10)
        break;
    end
end

% hold off

%儲存模型在子目錄model中，注：這裡目錄下必須有model資料夾不然報錯mkdir('model');可以解決
file_name=sprintf('model\\w_adagrad_n%d_lr%d.mat',i,lr);%檔名標準化;
save(file_name,'w');%儲存w向量儲存模型

clear detal_y  i iteration lr sum_grad

結果比較和處理

讀取test.csv

BTW:同樣的這裡NR全部已經改為0了

%讀取testing data ;test.csv中有240天的資料
x_test=ones(240,18*9+1);
for i=0:239
    for j=0:8%讀取前9小時資料 資料在CSV檔案的2到10列
        x_test(i+1,j*18+1:j*18+18)=(csvread('test.csv',i*18,j+2,[i*18,j+2, i*18+17,j+2 ]))';
    end
end
clear i j

結果儲存以及與實際值的比較

%讀取模型資料
w=open(file_name);%讀取的是struct結構
w=w.w;%將結構中的數匯入w中

%資料估計
y_test=x_test*w;

%儲存資料
y_file_name=strrep(file_name,'w','y'); %標準化輸出檔名
save(y_file_name,'y_test');

%計算估計值和實際值得方差
%讀取真實資料
y_real=csvread('ans.csv',1,1,[1,1,240,1]);
variance=sum( (y_real-y_test).^2 )/240;%計算方差

% variance_file_name=strrep(y_file_name,'y','variance'); %標準化輸出檔名
% variance_file_name=strrep(variance_file_name,'mat','txt');%標準化輸出檔名
fid=fopen('model\variance.csv','a+');
data_name=strrep(y_file_name,'y','variance');%標準化輸出
data_name=strrep(data_name,'.mat',',');%標準化輸出
data_name=strrep(data_name,'model\','\r\n');%標準化輸出
fprintf(fid,data_name);
fprintf(fid,'%.3f',variance);
fclose(fid);

clear y_file_name data_name fid