基于LSTM的一維數(shù)據(jù)擬合擴(kuò)展

一、引(fei)言(hua)

我在做Sri Lanka生態(tài)系統(tǒng)服務(wù)價值計算時，中間遇到了一點(diǎn)小問題。從世界糧農(nóng)組織(FAO)上獲得Sri Lanka主要農(nóng)作物產(chǎn)量和價格數(shù)據(jù)時，其中的主要作物Sorghum僅有2001-2006年的數(shù)據(jù)，而Millet只有2001-2005,2020-2021這樣的間斷數(shù)據(jù)。雖然說可以直接剔除這種過分缺失的數(shù)據(jù)，但這無疑會對生態(tài)因子的計算造成重大影響。所以我想要不要整個函數(shù)把他擬合一下，剛好Maize和Rice有2001-2021的完備數(shù)據(jù)，于是，這個文檔就這樣誕生了。

二、數(shù)據(jù)

數(shù)據(jù)來自FAO，考慮到可能有同學(xué)想要跟著嘗試一下，這里給出用到的數(shù)據(jù)。

作物產(chǎn)量

作物價格

2.1 數(shù)據(jù)探查

我們讀取數(shù)據(jù)，并進(jìn)行簡單的統(tǒng)計量查看。如果要進(jìn)一步深入研究數(shù)據(jù)分布及可視化，可以看看我的這篇文章

import pandas as pd

path=r"YourPath"

yield_=pd.read_csv(path+r"\yield.csv")
pp_=pd.read_csv(path+r"\Producer Prices.csv")

yield_.head()

需要用到的屬性只有Item,Year,Unit,Value

所以我們做這樣的處理：

yield_=yield_[["Item","Year","Unit","Value"]]

可以看到有些數(shù)據(jù)是從1961年開始的，太舊了就不用了，我們從2001年開始。

yield_=yield_[yield_["Year"]>2000]

同樣，我們來看看pp_的情況：

pp_.head()

pp_=pp_[["Item","Year","Value","Element"]]
pp_=pp_[pp_["Year"]>2000]

實際上，在這個數(shù)據(jù)里，產(chǎn)量已經(jīng)沒有問題了。我們只需要做一個簡單的處理：

yield_.groupby("Item").mean()["Value"]/10 #轉(zhuǎn)為千克

便可拿到每種作物近二十年的平均產(chǎn)量。

好了現(xiàn)在大問題出現(xiàn)在價值上，我們從下往上看就知道了：

pp_.tail(10)

高粱只有2006年的，那有沒有辦法利用現(xiàn)成的數(shù)據(jù)將其擴(kuò)展呢？

實際上，這類擬合問題有很多種解決方案，但是本問題涉及到時間，之前時間段的因子，以及可能的周期性，都會增加擬合的復(fù)雜性。所以，在這里我們采用LSTM來填充數(shù)據(jù)。

三、模型構(gòu)建

在本小節(jié)，我們將比較傳統(tǒng)一維CNN與RNN在結(jié)果上的異同。

一般做一維RNN時，可以指定一個時間窗口，比如用2006,2007,2008年的數(shù)據(jù)，推理2009年的數(shù)據(jù)，用2007,2008,2009年推理2010年。

我們現(xiàn)在要用之前處理好的pp_c數(shù)據(jù)中的玉米產(chǎn)量，來預(yù)測高粱產(chǎn)量。所以第一步就是將其轉(zhuǎn)化為torch接受的格式。

別忘記導(dǎo)入模塊：

import torch
import torch.nn as nn
from torch.nn import functional as F

x=pp_c[pp_c['Item']=="Maize (corn)"]['Value']
x=torch.FloatTensor(x)

之前寫數(shù)據(jù)迭代器的時候，除了可以繼承自torch.utils.data.DataLoader，也可以是任意的可迭代對象。這里我們可以簡單的設(shè)置一個類：

# 設(shè)置迭代器
class MyDataSet(object):
    def __init__(self,seq,ws=6):
        # ws是滑動窗口大小
        self.ori=[i for i in seq[:ws]]
        self.label=[i for i in seq[ws:]]
        self.reset()
        self.ws=ws

    def set(self,dpi):
        # 添加數(shù)據(jù)
        self.x.append(dpi)
        
    def reset(self):
        # 初始化
        self.x=self.ori[:]
        
    def get(self,idx):
        return self.x[idx:idx+self.ws],self.label[idx]
    
    def __len__(self):
        return len(self.x)

哦這邊提一下，有兩種方式，一種是用原始數(shù)據(jù)做預(yù)測，一種是用預(yù)測數(shù)據(jù)做預(yù)測，可能有點(diǎn)抽象，下面舉個例子。

假設(shè)$A=[a1,a2,a3,a4,a5,a6]$，時間窗口大小為3。

用原始數(shù)據(jù)做預(yù)測，那么輸入值為：$a1,a2,a3$，得到的結(jié)果將與$a4$做比較。下一輪輸入為$a2,a3,a4$，得到的結(jié)果將與$a5$做比較。

而用預(yù)測的數(shù)據(jù)做預(yù)測，第一輪輸入值為$a1,a2,a3$，得到的結(jié)果是$b4$，在與$a4$做比較后，下一輪的輸入為$a2,a3,b4$，會出現(xiàn)如下情況：

輸入數(shù)據(jù)為$b4,b5,b6$。

我們現(xiàn)在舉的例子是用預(yù)測的數(shù)據(jù)做預(yù)測。當(dāng)然，最后也會給出一個用原始數(shù)據(jù)做預(yù)測的版本，那個版本相對簡單。

ws=6 # 全局時間窗口
train_data=MyDataSet(x,ws)

網(wǎng)絡(luò)的架構(gòu)如下：

   
class Net3(nn.Module):
    def __init__(self,in_features=54,n_hidden1=128,n_hidden2=256,n_hidden3=512,out_features=7):
        super(Net3, self).__init__()
        self.flatten=nn.Flatten()
        self.hidden1=nn.Sequential(
            nn.Linear(in_features,n_hidden1,False),
           
            nn.ReLU()
        )
        self.hidden2=nn.Sequential(
            nn.Linear(n_hidden1,n_hidden2),

            nn.ReLU()
        )
        self.hidden3=nn.Sequential(
            nn.Linear(n_hidden2,n_hidden3),

            nn.ReLU()
        )
        self.out=nn.Sequential(nn.Linear(n_hidden3,out_features))

    def forward(self,x):
        x=self.flatten(x)
        x=self.hidden2(self.hidden1(x))
        x=self.hidden3(x)
        return self.out(x)



class CNN(nn.Module):
    def __init__(self, output_dim=1,ws=6):
        super(CNN, self).__init__()
        self.relu = nn.ReLU(inplace=True)
        self.conv1 = nn.Conv1d(ws, 64, 1)
        self.lr = nn.LeakyReLU(inplace=True)
        self.conv2 = nn.Conv1d(64, 128, 1)

        self.bn1, self.bn2 = nn.BatchNorm1d(64), nn.BatchNorm1d(128)
        self.bn3, self.bn4 = nn.BatchNorm1d(1024), nn.BatchNorm1d(128)
        self.flatten = nn.Flatten()
        self.lstm1 = nn.LSTM(128, 1024)
        self.lstm2 = nn.LSTM(1024, 256)
        self.lstm3=nn.LSTM(256,512)
        self.fc = nn.Linear(512, 512)
        self.fc4=nn.Linear(512,256)
        self.fc1 = nn.Linear(256, 64)
        self.fc3 = nn.Linear(64, output_dim)

    @staticmethod
    def reS(x):
        return x.reshape(-1, x.shape[-1], x.shape[-2])

    def forward(self, x):
        x = self.reS(x)
        x = self.conv1(x) 
        x = self.lr(x)

        x = self.conv2(x) 
        x = self.lr(x)

        x = self.flatten(x)

        # LSTM部分
        x, h = self.lstm1(x)
        x, h = self.lstm2(x)
        x,h=self.lstm3(x)
        x, _ = h

        x = self.fc(x.reshape(-1, ))
        x = self.relu(x)
        x = self.fc4(x)
        x = self.relu(x)
        x = self.fc1(x)
        x = self.relu(x)
        x = self.fc3(x)

        return x

Net3主要是一維卷積，CNN加入了LSTM結(jié)構(gòu)。至于名字，是隨便取的…跟內(nèi)容并無關(guān)系。

def Train(model,train_data,seed=1):
    device="cuda" if torch.cuda.is_available() else "cpu"
    model=model.to(device)
    Mloss=100000
    path=r"YourPath\%s.pth"%seed
    # 設(shè)置損失函數(shù),這里使用的是均方誤差損失
    criterion = nn.MSELoss()
    # 設(shè)置優(yōu)化函數(shù)和學(xué)習(xí)率lr
    optimizer=torch.optim.Adam(model.parameters(),lr=1e-5,betas=(0.9,0.99),
                               eps=1e-07,weight_decay=0)
    # 設(shè)置訓(xùn)練周期
    epochs =3000
    criterion=criterion.to(device)
    model.train()
    for epoch in range(epochs):
        total_loss=0

        for i in range(len(x)-ws):
            # 每次更新參數(shù)前都梯度歸零和初始化
            seq,y_train=train_data.get(i) # 從我們的數(shù)據(jù)集中拿出數(shù)據(jù)
            seq,y_train=torch.FloatTensor(seq),torch.FloatTensor([y_train])
            seq=seq.unsqueeze(dim=0)
            seq,y_train=seq.to(device),y_train.to(device)

            optimizer.zero_grad()
            # 注意這里要對樣本進(jìn)行reshape，
            # 轉(zhuǎn)換成conv1d的input size（batch size, channel, series length）
            y_pred = model(seq)
            loss = criterion(y_pred, y_train)
            loss.backward()
            train_data.set(y_pred.to("cpu").item()) # 再放入預(yù)測數(shù)據(jù)
            optimizer.step()
            total_loss+=loss

        train_data.reset()
        if total_loss.tolist()<Mloss:
            Mloss=total_loss.tolist()
            torch.save(model.state_dict(),path)
            print("Saving")
        print(f'Epoch: {epoch+1:2} Mean Loss: {total_loss.tolist()/len(train_data):10.8f}')
    return model

正常訓(xùn)練就OK

d=CNN(ws=ws)
Train(d,train_data,4)

平均損失在10點(diǎn)左右，還有很大優(yōu)化空間。當(dāng)然我們這里只是舉個非常簡單的例子，就是個baseline

checkpoint=torch.load(r"YourPath\4.pth")
d.load_state_dict(checkpoint) # 加載最佳參數(shù)
d.to("cpu")

四、結(jié)果可視化

我們這里用到Pyechart進(jìn)行可視化。

from pyecharts.charts import *
from pyecharts import options as opts
from pyecharts.globals import CurrentConfig

pre,ppre=[i.item() for i in x[:ws]],[]
# pre 是用原始數(shù)據(jù)做預(yù)測
# ppre 用預(yù)測數(shù)據(jù)做預(yù)測
for i in range(len(x)-ws+1):
    ppre.append(d(torch.FloatTensor(x[i:i+ws]).unsqueeze(dim=0)))
    pre.append(d(torch.FloatTensor(pre[-ws:]).unsqueeze(dim=0)).item())

l=Line()
l.add_xaxis([i for i in range(len(x))])
l.add_yaxis("Original Data",x.tolist())
l.add_yaxis("Pred Data(Using Raw Datas)",x[:ws].tolist()+[i.item() for i in ppre])
l.add_yaxis("Pred Data(Using Pred Datas)",pre)
l.set_series_opts(label_opts=opts.LabelOpts(is_show=False))
l.set_global_opts(title_opts=opts.TitleOpts(title='LSTM CNN'))

l.render_notebook()

根據(jù)時間窗口的不同，可以得到不同的結(jié)果。

ws=4

ws=5

ws=6

從結(jié)果上來看，時間窗口越大越好。但是這里我們只能到六了，再大就不禮貌了。(高粱只有六個節(jié)點(diǎn)的數(shù)據(jù))。

至于驗證，我們可以選Rice做驗證：

x=torch.FloatTensor(pp_c[pp_c['Item']=="Rice"]['Value'].tolist())
pre,ppre=[i.item() for i in x[:ws]],[]
for i in range(len(x)-ws+1):
    ppre.append(d(torch.FloatTensor(x[i:i+ws]).unsqueeze(dim=0)))
    pre.append(d(torch.FloatTensor(pre[-ws:]).unsqueeze(dim=0)).item())
l=Line()
l.add_xaxis([i for i in range(len(x))])
l.add_yaxis("Original Data",x.tolist())
l.add_yaxis("Pred Data(Using Raw Datas)",x[:ws].tolist()+[i.item() for i in ppre])
l.add_yaxis("Pred Data(Using Pred Datas)",pre)
l.set_series_opts(label_opts=opts.LabelOpts(is_show=False))
l.set_global_opts(title_opts=opts.TitleOpts(title='LSTM CNN'))

l.render_notebook()

可以發(fā)現(xiàn)，用預(yù)測做預(yù)測的結(jié)果，基本上不會差太多，那也就意味著，我們可以對高粱進(jìn)行預(yù)測啦！不過在這之前，我們可以看看用原始數(shù)據(jù)做訓(xùn)練的結(jié)果：

時間窗口一樣為6，可以看到在黑線貼合的非常好，但是面對大量缺失的數(shù)據(jù)，精度就遠(yuǎn)不如用預(yù)測數(shù)據(jù)做預(yù)測的結(jié)果了。

此外，這是用CNN做的結(jié)果

我們可以發(fā)現(xiàn)LSTM的波動要比CNN好，CNN后面死水一潭，應(yīng)該是梯度消失導(dǎo)致的，前面信息沒有了，后面信息又是自個構(gòu)造的，這就導(dǎo)致了到后面變成了線性情況。

那么最后的最后，就是預(yù)測高粱產(chǎn)量了：

pre_data=pp_c[pp_c['Item']=='Sorghum']['Value'].tolist()
l=pre_data[:]
for i in range(len(x)-ws+1):
    l.append(d(torch.FloatTensor(l[-ws:]).unsqueeze(dim=0)).item())
L=Line()
L.add_xaxis([i for i in range(len(x))])
L.add_yaxis("Pred",l)
L.set_series_opts(label_opts=opts.LabelOpts(is_show=False))
L.set_global_opts(title_opts=opts.TitleOpts(title='sorghum production forecasts')
                            
                             )

L.render_notebook()
l.to_csv("path")

文章來源地址http://www.zghlxwxcb.cn/news/detail-658062.html

到了這里，關(guān)于【項目實踐】基于LSTM的一維數(shù)據(jù)擴(kuò)展與預(yù)測的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！