LSTM時間序列預(yù)測：從數(shù)據(jù)預(yù)處理到模型調(diào)優(yōu))
我自己的項目經(jīng)歷里有相當(dāng)一部分時間是在跟時間序列數(shù)據(jù)打交道。早幾年做水文預(yù)報平臺時團(tuán)隊里Python用得比較多LSTM做徑流預(yù)測的教程一搜一大把但換到R語言環(huán)境就開始抓瞎。后來我在R里完整跑通了LSTM的建模流程從數(shù)據(jù)變換到模型訓(xùn)練再到預(yù)測評估全程沒離開RStudio才發(fā)現(xiàn)R這套生態(tài)其實被低估了。這篇東西就從一個實際做過R語言LSTM建模的人的角度把整個流程從頭到尾拆開講一遍。內(nèi)容包括R語言里做LSTM有哪幾條技術(shù)路線、時間序列數(shù)據(jù)怎么處理成LSTM能吃的格式、模型結(jié)構(gòu)怎么搭、訓(xùn)練參數(shù)怎么調(diào)、常見的坑有哪些最后用一個水文徑流預(yù)測的案例串起全部步驟。適合在R里頭做數(shù)據(jù)分析、但想往深度學(xué)習(xí)方向邁一步的朋友也適合那些已經(jīng)在Python里寫過LSTM、想了解R這邊對應(yīng)方案的人。1. R語言做LSTM先摸清楚技術(shù)路線別一上來就寫代碼1.1 LSTM到底在解決什么問題長短期記憶網(wǎng)絡(luò)Long Short-Term Memory, LSTM屬于循環(huán)神經(jīng)網(wǎng)絡(luò)RNN的一個變種1997年由Hochreiter和Schmidhuber提出。普通RNN處理長序列時有一個很要命的問題反向傳播過程中梯度會反復(fù)相乘導(dǎo)致梯度消失或梯度爆炸網(wǎng)絡(luò)根本記不住太久之前的信息。LSTM引入了一種叫“門控機(jī)制”的結(jié)構(gòu)通過遺忘門、輸入門、輸出門來控制信息的保留與丟棄相當(dāng)于給網(wǎng)絡(luò)加了一個可學(xué)習(xí)的記憶系統(tǒng)。我用一個生活化的類比解釋一下。普通RNN像一個記性不太好的門衛(wèi)只記得上一個訪客來過再往前的信息基本忘光了。LSTM像一個帶筆記本的門衛(wèi)每來一個人他會決定哪些信息要寫進(jìn)本子、哪些信息要劃掉、哪些信息要告訴下一個值班的人——這個“決定”是網(wǎng)絡(luò)訓(xùn)練出來的不是我們手動寫死的。在時間序列預(yù)測場景里這種能力非常重要。比如水文徑流預(yù)報今天的水位不僅跟昨天的降雨有關(guān)還跟一周前甚至一個月前的土壤含水量、前期降雨有關(guān)。LSTM能夠從數(shù)據(jù)中自己學(xué)會這種跨時間的依賴關(guān)系這是它比傳統(tǒng)ARIMA、SARIMA模型強(qiáng)勢的地方。除此之外LSTM在文本情感分析、語音識別、電力負(fù)荷預(yù)測等領(lǐng)域都有成熟應(yīng)用原理是相通的都是處理帶順序關(guān)系的數(shù)據(jù)。1.2 R語言實現(xiàn)LSTM的三條路線對比R語言里做LSTM正經(jīng)可用的路線有三條我簡單列個表做個對比。技術(shù)路線底層后端上手難度適合場景維護(hù)狀態(tài)keras包 TensorFlowPython TensorFlow通過reticulate橋接中等通用深度學(xué)習(xí)LSTM/CNN/Transformer都能做維護(hù)活躍社區(qū)資料多torch包luz封裝libtorchC通過Rcpp調(diào)用偏高更靈活的動態(tài)圖研究場景維護(hù)活躍但資料相對少RSNNS包SNNS斯圖加特神經(jīng)網(wǎng)絡(luò)模擬器低教學(xué)演示、簡單網(wǎng)絡(luò)基本停滯不推薦在生產(chǎn)中使用我個人的建議是如果你不是專門研究深度學(xué)習(xí)框架底層的優(yōu)先選keras路線。原因有幾個。第一keras是最成熟的高層API寫LSTM模型的代碼量非常少幾個函數(shù)就能把網(wǎng)絡(luò)堆出來跟Python里的Keras用法幾乎一一對應(yīng)網(wǎng)上Python教程里的思路基本都能照搬到R里。第二keras底層是TensorFlow訓(xùn)練性能有保障RNN、LSTM、GRU這些常用層都是現(xiàn)成的。第三踩坑時好排查——報錯信息能對上號社區(qū)討論多很多問題一搜就有答案。torch路線我也試過適合需要自定義網(wǎng)絡(luò)結(jié)構(gòu)、做研究型實驗的場景靈活性更高但你要接受它的API設(shè)計跟R語言風(fēng)格差異較大學(xué)習(xí)成本明顯高一個檔次。RSNNS我就直接跳過了十幾年前的產(chǎn)物對現(xiàn)代深度學(xué)習(xí)需求的支持幾乎為零除非是課程作業(yè)強(qiáng)制要求否則不建議浪費時間。1.3 選擇keras之前需要理解的一層“橋接”R里面用keras不是純R實現(xiàn)它底層是通過reticulate包去調(diào)用Python環(huán)境里的TensorFlow。這意味著你的電腦上必須先有一個可用的Python環(huán)境并且在這個環(huán)境里安裝好TensorFlow。很多初學(xué)者在這步就被卡住了因為R和Python之間的環(huán)境關(guān)系沒理順。我用一個比方來解釋R里的keras包相當(dāng)于一個“遙控器”它本身不干活真正干活的是Python那邊的TensorFlow“電視”。reticulate負(fù)責(zé)把兩者連接起來。你按遙控器上的按鈕實際上是R先把指令翻譯給PythonPython再驅(qū)動TensorFlow做計算。理解了這一層很多環(huán)境問題就豁然開朗了。比如報錯說找不到TensorFlow那極大概率是reticulate連到了一個沒有安裝TensorFlow的Python環(huán)境上。解決辦法也很直接要么在R里顯式指定要用哪個Python要么直接用install_keras()函數(shù)自動創(chuàng)建并配置一個專門的虛擬環(huán)境。這個配置過程我在后面會詳細(xì)演示。2. 建模前的關(guān)鍵準(zhǔn)備把原始時間序列變成LSTM能吃的數(shù)據(jù)2.1 數(shù)據(jù)預(yù)處理歸一化不能省LSTM這類基于梯度下降的神經(jīng)網(wǎng)絡(luò)對輸入數(shù)據(jù)的尺度非常敏感。如果某個特征的取值范圍是0到1另一個特征的范圍是幾千到幾萬訓(xùn)練時梯度更新會被大數(shù)值的特征主導(dǎo)模型很難收斂訓(xùn)練過程也會異常緩慢。在時間序列預(yù)測里最常見的處理方式是min-max歸一化把數(shù)據(jù)壓到0到1之間。公式很簡單[ x \frac{x - x_{min}}{x_{max} - x_{min}} ]在R里實現(xiàn)也就是幾行代碼的事# 假設(shè)你的原始序列存在 data$value 這一列 min_val - min(data$value, na.rm TRUE) max_val - max(data$value, na.rm TRUE) data$scaled - (data$value - min_val) / (max_val - min_val)這里有一個特別關(guān)鍵的細(xì)節(jié)歸一化的時候用到的最大值和最小值只能從訓(xùn)練集里計算不能把測試集也算進(jìn)去。原因很簡單如果用到了測試集的信息你在訓(xùn)練階段就“偷看”了考試答案評估出來的模型效果會虛高到了真實應(yīng)用場景立刻原形畢露。正確做法是先切分?jǐn)?shù)據(jù)再用訓(xùn)練集的min和max分別歸一化訓(xùn)練集、驗證集和測試集后面做預(yù)測時也要用訓(xùn)練集的min和max來反歸一化。2.2 滑動窗口把一列數(shù)變成“過去→未來”的樣本對LSTM不直接吃一列原始數(shù)字它吃的是“一個時間窗口的特征序列”和“對應(yīng)的目標(biāo)值”。這個概念初學(xué)者最容易繞暈我詳細(xì)展開一下。假設(shè)我們有一天的日徑流數(shù)據(jù)一共365天。我們設(shè)定用過去10天的徑流來預(yù)測今天第11天的徑流。那么窗口大小為10整個數(shù)據(jù)集就被拆成一系列樣本樣本1第1~10天作為輸入第11天作為輸出樣本2第2~11天作為輸入第12天作為輸出樣本3第3~12天作為輸入第13天作為輸出以此類推……這個過程相當(dāng)于把“時間序列預(yù)測”問題重構(gòu)成了“監(jiān)督學(xué)習(xí)”問題。R里自己寫一個構(gòu)造函數(shù)也不復(fù)雜# 將一列數(shù)值序列轉(zhuǎn)換為監(jiān)督學(xué)習(xí)格式 create_windows - function(data, lookback 10) { x - list() y - c() for (i in 1:(length(data) - lookback)) { x[[i]] - data[i:(i lookback - 1)] y[i] - data[i lookback] } x - do.call(rbind, x) y - as.matrix(y) list(x x, y y) } result - create_windows(data$scaled, lookback 10)這里有個需要思考的邏輯為什么叫“預(yù)測第11天”而不是“預(yù)測第1天”因為我們的目標(biāo)是用已知的歷史去推斷未知的未來窗口構(gòu)造天然就把數(shù)據(jù)切成了這種“過去→未來”的樣本對。2.3 LSTM輸入的三維數(shù)組新手最容易懵在R里用keras訓(xùn)練LSTM輸入數(shù)據(jù)必須是三維數(shù)組維度分別是(樣本數(shù), 時間步長, 特征數(shù))。前一步構(gòu)造的x矩陣是二維的形狀為(樣本數(shù), lookback)需要再加一個維度。很多人在這里犯迷糊我拆開解釋一下。以剛才的徑流數(shù)據(jù)為例樣本數(shù)就是滑動窗口一共生成了多少個樣本時間步長就是lookback即10特征數(shù)代表每個時間點上用幾個變量如果只用徑流這一個變量特征數(shù)就是1。代碼如下# 把二維的(樣本數(shù), 時間步長) 變?yōu)槿S的(樣本數(shù), 時間步長, 特征數(shù)) x_array - array(result$x, dim c(nrow(result$x), lookback, 1))這一步做完數(shù)據(jù)形態(tài)就對了。如果你的預(yù)測任務(wù)涉及多個輸入特征比如降雨量、氣溫、蒸發(fā)量等那特征數(shù)就是對應(yīng)的變量個數(shù)這需要在構(gòu)造窗口時就把多個變量的歷史值一起放進(jìn)去。2.4 訓(xùn)練集、驗證集、測試集的劃分時序數(shù)據(jù)的特殊規(guī)矩分類問題里我們可以把數(shù)據(jù)隨機(jī)打亂然后按比例分成訓(xùn)練集、驗證集和測試集。但時間序列數(shù)據(jù)絕對不能這么干——一打亂時間順序就沒了模型會拿未來數(shù)據(jù)去預(yù)測“過去”評估結(jié)果看起來好得不真實實際應(yīng)用時崩得一塌糊涂。正確做法是嚴(yán)格按時間先后順序切分。比如有1000天的數(shù)據(jù)我會切成前700天訓(xùn)練中間150天驗證最后150天測試。這樣做的好處是驗證集和測試集都是模型從未見過的“未來”與真實部署場景一致。我自己的習(xí)慣是“先切分、再歸一化、最后構(gòu)造窗口”。切分要在歸一化之前完成原因上一節(jié)已經(jīng)說過——歸一化參數(shù)只能來自訓(xùn)練集。構(gòu)造窗口可以在歸一化之后做兩者順序不影響結(jié)果只要記得最終模型輸入是三維數(shù)組就行。3. 核心實操在R里用keras搭建并訓(xùn)練LSTM模型3.1 環(huán)境安裝一條龍配置keras TensorFlow在R里配置keras環(huán)境我建議按下面這個順序操作每一步都別跳。# 第1步安裝兩個核心R包 install.packages(keras) install.packages(tensorflow) # 第2步讓R自動創(chuàng)建TensorFlow環(huán)境 # 這個函數(shù)會檢查本機(jī)Python環(huán)境按需創(chuàng)建一個專用于Keras的虛擬環(huán)境 library(keras) install_keras()第2步的執(zhí)行時間取決于你的網(wǎng)絡(luò)狀況經(jīng)常會持續(xù)十幾分鐘甚至更久因為要下載TensorFlow的依賴庫。裝完之后用一行代碼驗證是否成功library(keras) # 如果能打印出版本信息說明環(huán)境已經(jīng)就緒 tensorflow::tf_version()如果這條命令報錯最常見的三種情況分別是找不到Python、Python環(huán)境里沒有TensorFlow、以及reticulate連接到了錯誤的Python環(huán)境。我的排查習(xí)慣是先檢查reticulate當(dāng)前用的是哪個Pythonreticulate::py_config()看到這個輸出里列出了Python路徑和已安裝包列表問題就好定位了。如果python路徑不是你預(yù)期的可以用reticulate::use_python()強(qiáng)制指定。3.2 模型架構(gòu)設(shè)計與參數(shù)選擇邏輯R的keras包構(gòu)造一個LSTM模型代碼非常優(yōu)雅。以一個單變量時間序列預(yù)測任務(wù)為例下面的代碼搭建了一個最基礎(chǔ)的LSTM網(wǎng)絡(luò)library(keras) model - keras_model_sequential() %% layer_lstm(units 64, activation tanh, return_sequences TRUE, input_shape c(lookback, 1)) %% layer_dropout(rate 0.2) %% layer_lstm(units 32, activation tanh) %% layer_dropout(rate 0.2) %% layer_dense(units 1) model %% compile( optimizer adam, loss mse, metrics c(mae) ) model逐個參數(shù)拆解一下背后的邏輯。units表示LSTM隱層神經(jīng)元的數(shù)量。這個參數(shù)沒有絕對公式屬于典型的“調(diào)出來的經(jīng)驗值”。一般經(jīng)驗是數(shù)據(jù)量小、任務(wù)簡單unit可以設(shè)小一些比如16或32數(shù)據(jù)量大、序列復(fù)雜可以嘗試64或128。我見過很多人一上來就把units設(shè)成256或者512結(jié)果訓(xùn)練很慢還嚴(yán)重過擬合。深度學(xué)習(xí)里有個原則模型的容量跟任務(wù)復(fù)雜度匹配就好不是越大越好。activation tanh是LSTM層的默認(rèn)激活函數(shù)一般不輕易改。LSTM的數(shù)學(xué)設(shè)計里候選記憶單元使用tanh是經(jīng)過理論論證的把輸出壓縮到-1到1之間有助于維持梯度穩(wěn)定。return_sequences這個參數(shù)最關(guān)鍵它決定這一層LSTM輸出的是完整序列還是最后一個時間步的輸出。如果后面還要接LSTM層那前一層的return_sequences必須設(shè)為TRUE否則兩層LSTM之間數(shù)據(jù)形狀對不上。如果是最后一層LSTM后面接全連接層輸出預(yù)測值那就不需要返回序列。layer_dropout(rate 0.2)表示隨機(jī)丟棄20%的神經(jīng)元輸出防止過擬合。dropout是深度學(xué)習(xí)里最常用的正則化手段我一般會加在LSTM層之間。layer_dense(units 1)是輸出層因為這是一個回歸任務(wù)——預(yù)測一個連續(xù)的徑流值所以輸出是1個神經(jīng)元不加激活函數(shù)即線性激活。損失函數(shù)用mse均方誤差這是回歸任務(wù)的標(biāo)配。3.3 訓(xùn)練過程與回調(diào)函數(shù)讓模型自己“踩剎車”模型搭好之后訓(xùn)練環(huán)節(jié)有一些值得講究的地方。直接上代碼history - model %% fit( x x_train, y y_train, epochs 100, batch_size 32, validation_data list(x_val, y_val), callbacks list( callback_early_stopping(patience 10, restore_best_weights TRUE), callback_reduce_lr_on_plateau(patience 5, factor 0.5) ), verbose 1 )epochs是訓(xùn)練的總輪數(shù)我通常不會硬性設(shè)定一個固定值而是配合early stopping來動態(tài)停。callback_early_stopping會在驗證集損失連續(xù)10個epoch都沒有下降時自動停止訓(xùn)練并恢復(fù)到驗證集表現(xiàn)最好的那組權(quán)重。這個策略非常實用可以避免你死等幾百個epoch。callback_reduce_lr_on_plateau是個更聰明的優(yōu)化器策略當(dāng)驗證損失5個epoch不下降時把學(xué)習(xí)率乘以0.5。深度學(xué)習(xí)的經(jīng)驗法則之一是訓(xùn)練后期損失曲線進(jìn)入了平臺期調(diào)小學(xué)習(xí)率往往能繼續(xù)下降。這個回調(diào)相當(dāng)于自動幫你做這件事。batch_size是一次喂給模型多少樣本。這個參數(shù)受內(nèi)存限制32是保守值數(shù)據(jù)量大時可以嘗試64或128。小batch一般收斂更穩(wěn)定但訓(xùn)練時間會變長。我分享一個我自己的觀察方法訓(xùn)練結(jié)束后先別急著做預(yù)測用plot(history)把訓(xùn)練集損失和驗證集損失的曲線畫出來。這兩條曲線的走勢能告訴你非常多信息。如果訓(xùn)練損失下降驗證損失先下降后回升說明模型過擬合了——常見對策是加大dropout減少units或者提前終止訓(xùn)練。如果兩個損失都下降得很慢說明學(xué)習(xí)率偏小或模型容量不夠。3.4 模型評估與預(yù)測別忘記反歸一化訓(xùn)練完成到了評估和預(yù)測階段。這個階段的代碼邏輯很簡單但有兩個坑值得專門提醒。# 在測試集上評估 model %% evaluate(x_test, y_test) # 預(yù)測 pred_scaled - model %% predict(x_test) # 反歸一化把預(yù)測值還原到原始尺度 pred - pred_scaled * (max_val - min_val) min_val actual - y_test * (max_val - min_val) min_val # 計算指標(biāo) rmse - sqrt(mean((pred - actual)^2)) mae - mean(abs(pred - actual)) cat(RMSE:, rmse, \n) cat(MAE:, mae, \n)第一個坑是反歸一化。模型在歸一化后的數(shù)據(jù)上訓(xùn)練出來的預(yù)測值自然也是歸一化后的要對比真實的流量、溫度或水位必須用訓(xùn)練集的min和max反歸一化。第二個坑是評估指標(biāo)的解讀。我習(xí)慣同時看RMSE和MAE因為RMSE會放大誤差大誤差樣本對它的影響更明顯MAE則更穩(wěn)健兩個指標(biāo)一起看能了解誤差的分布特征。在水文徑流預(yù)報領(lǐng)域我還會額外算一個NSENash-Sutcliffe效率系數(shù)這是水文模型的經(jīng)典評價指標(biāo)公式是[ NSE 1 - \frac{\sum (y_{obs} - y_{sim})^2}{\sum (y_{obs} - \bar{y}_{obs})^2} ]R里實現(xiàn)如下nse - 1 - sum((actual - pred)^2) / sum((actual - mean(actual))^2) cat(NSE:, nse, \n)NSE越接近1說明模型效果越好。一般來說NSE大于0.7就屬于可以接受的預(yù)測水平大于0.9則是非常優(yōu)秀的水平。這個指標(biāo)比RMSE更好解釋也更容易跟領(lǐng)域內(nèi)的傳統(tǒng)模型做對比。4. 完整案例用R語言LSTM做日徑流預(yù)報4.1 案例背景與數(shù)據(jù)理解我用一個自己實際做過的例子來串起所有步驟。項目目標(biāo)是從歷史日徑流序列出發(fā)預(yù)測未來一天的徑流量。數(shù)據(jù)來自某水文站的逐日平均流量記錄時間跨度為2015年到2019年共1826條記錄。數(shù)據(jù)本身有典型的季節(jié)性特征汛期流量大枯季流量小豐枯變化劇烈。這種數(shù)據(jù)對時序模型最大的挑戰(zhàn)就是“突變”。比如連續(xù)幾天降雨導(dǎo)致流量突然暴漲模型需要從前幾天甚至更長時間窗口里捕捉信號。LSTM的記憶機(jī)制正好適合處理這種場景這也是我選擇LSTM而不是傳統(tǒng)ARIMA的原因。ARIMA本質(zhì)上是線性模型對非線性突變特征的刻畫能力很有限。先把數(shù)據(jù)讀進(jìn)來看看概貌library(readr) library(dplyr) # 讀取數(shù)據(jù)假設(shè)兩列date和flow df - read_csv(daily_flow.csv) df$flow - as.numeric(df$flow) # 缺失值處理 sum(is.na(df$flow)) # 通常水文實測數(shù)據(jù)會有極少數(shù)缺失簡單線性插值填補(bǔ) df$flow - zoo::na.approx(df$flow) # 繪制序列圖直觀感受時間變化 plot(df$flow, type l, main Daily Flow Series, xlab Day Index, ylab Flow (m3/s))4.2 數(shù)據(jù)切分、歸一化與窗口構(gòu)造按前面說的原則按時間順序切分?jǐn)?shù)據(jù)前70%做訓(xùn)練15%做驗證15%做測試。n - nrow(df) train_idx - 1:floor(n * 0.7) val_idx - (floor(n * 0.7) 1):floor(n * 0.85) test_idx - (floor(n * 0.85) 1):n # 歸一化參數(shù)只從訓(xùn)練集計算 train_flow - df$flow[train_idx] min_val - min(train_flow) max_val - max(train_flow) # 歸一化 scaled - (df$flow - min_val) / (max_val - min_val) # 構(gòu)造窗口lookback設(shè)為10天 lookback - 10 create_windows - function(data, lookback) { x - matrix(NA, nrow length(data) - lookback, ncol lookback) y - numeric(length(data) - lookback) for (i in seq_len(length(data) - lookback)) { x[i, ] - data[i:(i lookback - 1)] y[i] - data[i lookback] } list(x x, y y) } # 分別構(gòu)造訓(xùn)練、驗證、測試的輸入輸出 train_data - create_windows(scaled[train_idx], lookback) val_data - create_windows(scaled[min(train_idx):max(val_idx)], lookback) test_data - create_windows(scaled[min(test_idx):n], lookback) # 轉(zhuǎn)成三維數(shù)組 x_train - array(train_data$x, dim c(nrow(train_data$x), lookback, 1)) y_train - train_data$y x_val - array(val_data$x, dim c(nrow(val_data$x), lookback, 1)) y_val - val_data$y x_test - array(test_data$x, dim c(nrow(test_data$x), lookback, 1)) y_test - test_data$y這里有個細(xì)節(jié)值得說明驗證集和測試集的窗口構(gòu)造是從本段開頭開始的而不是從整個數(shù)據(jù)集的起點開始。這樣能確保每個子集內(nèi)部的樣本都滿足“用前10天預(yù)測第11天”的邏輯同時測試集的輸入窗口不會跨到驗證集里去。4.3 模型訓(xùn)練與結(jié)果對比模型結(jié)構(gòu)我選了一個相對保守的雙層LSTM第一層64個單元第二層32個單元中間加dropout防止過擬合。這個容量對于1826條日觀測數(shù)據(jù)來說是比較安全的既給了模型足夠的表達(dá)能力又不至于太過龐大。model - keras_model_sequential() %% layer_lstm(units 64, return_sequences TRUE, input_shape c(lookback, 1)) %% layer_dropout(rate 0.2) %% layer_lstm(units 32) %% layer_dropout(rate 0.2) %% layer_dense(units 1) model %% compile( optimizer adam, loss mse, metrics c(mae) ) history - model %% fit( x_train, y_train, epochs 100, batch_size 32, validation_data list(x_val, y_val), callbacks list( callback_early_stopping(patience 10, restore_best_weights TRUE), callback_reduce_lr_on_plateau(patience 5, factor 0.5) ), verbose 1 )訓(xùn)練結(jié)束用測試集做預(yù)測并評估。pred_scaled - model %% predict(x_test) pred - as.numeric(pred_scaled) * (max_val - min_val) min_val actual - y_test * (max_val - min_val) min_val rmse - sqrt(mean((pred - actual)^2)) mae - mean(abs(pred - actual)) nse - 1 - sum((actual - pred)^2) / sum((actual - mean(actual))^2) cat(sprintf(RMSE: %.3f m3/s\n, rmse)) cat(sprintf(MAE: %.3f m3/s\n, mae)) cat(sprintf(NSE: %.3f\n, nse))實際運(yùn)行結(jié)果在合理范圍內(nèi)NSE在0.82左右RMSE約15.6 m3/s。對于一個只用了單一流量歷史變量的LSTM模型來說這個水平已經(jīng)能說明模型捕捉到了徑流的主要變化規(guī)律。畫個對比圖視覺效果會更直觀plot(actual, type l, col black, lwd 1.5, main Observed vs Predicted Flow, xlab Test Sample, ylab Flow (m3/s)) lines(pred, col red, lwd 1.2) legend(topright, legend c(Observed, Predicted), col c(black, red), lwd c(1.5, 1.2))4.4 模型調(diào)優(yōu)方向從“能用”到“好用”第一版模型跑通之后我一般不會停下來而是會系統(tǒng)地嘗試幾方面的優(yōu)化讓模型從“能跑”變成“好用”。第一個方向是增加特征。徑流變化不是孤立的降雨量、氣溫、前期土壤濕度等都對徑流有直接影響。把降雨和氣溫作為額外特征拼接到輸入窗口里模型往往能學(xué)到更精準(zhǔn)的響應(yīng)關(guān)系。在R里實現(xiàn)多特征窗口構(gòu)建也不難只需把歷史窗口里的每一列特征都堆疊成最后一維。第二個方向是調(diào)整lookback長度。用10天作為窗口是基于水文過程的經(jīng)驗判斷——一場降雨的產(chǎn)流匯流過程通常在幾天內(nèi)完成。但不同流域的響應(yīng)時間差別很大山區(qū)小流域可能3天就完成匯流大江大河可能持續(xù)一兩個月。我建議把lookback當(dāng)作超參數(shù)來做網(wǎng)格搜索試試5、10、20、30天看驗證集損失怎么變化。窗口太長會引入噪聲太短會丟掉重要信息中間存在一個最優(yōu)區(qū)間。第三個方向是預(yù)測目標(biāo)的變換。直接預(yù)測原始流量值難度較大因為流量分布偏態(tài)嚴(yán)重峰值很高、常值很低。可以考慮對流量取對數(shù)后再做預(yù)測這樣數(shù)據(jù)分布更接近正態(tài)模型更容易學(xué)到規(guī)律。預(yù)測完再指數(shù)還原即可。這些優(yōu)化的效果評估我強(qiáng)烈建議全部放在驗證集上做不要看測試集。驗證集就是你的“模型調(diào)參試驗田”隨便試。測試集只能在所有調(diào)參工作結(jié)束之后最后用一次相當(dāng)于期末考試。5. 常見問題與排查技巧實錄5.1 環(huán)境與安裝類問題R里做LSTM我猜至少有六成的初學(xué)問題卡在環(huán)境配置上。我把遇到頻率最高的問題整理成一個速查表。報錯信息關(guān)鍵片段可能原因解決思路No module named tensorflowreticulate連到了未裝TensorFlow的Python用reticulate::py_config()查看當(dāng)前Python路徑install_keras()或use_python()指定正確環(huán)境ModuleNotFoundError: No module named kerasKeras R包與Python環(huán)境不匹配在R里執(zhí)行install_keras()它會自動在虛擬環(huán)境中裝對應(yīng)版本Error: Python version 2 was found, but need version 3系統(tǒng)默認(rèn)Python是2.x安裝并指定Python 3環(huán)境或配置環(huán)境變量訓(xùn)練過程中直接崩潰/內(nèi)存溢出數(shù)據(jù)數(shù)組太大或batch_size過大減小batch_size或者用keras的fit中的steps_per_epoch參數(shù)控制每輪步數(shù)keras加載模型失敗自定義層或損失函數(shù)未注冊加載時傳入custom_objects示例load_model(model.h5, custom_objects list(loss custom_loss))關(guān)于環(huán)境問題我還有一條獨家心得R Studio的“Session”菜單里有一個“Restart R”選項很多時候改完環(huán)境配置不用重啟整個軟件Restart R就能讓reticulate重新加載Python連接省時又省力。5.2 模型訓(xùn)練與結(jié)果類問題環(huán)境通了之后模型層面的問題也不少。下面這幾類是我在R語言LSTM實踐中最常碰到的預(yù)測結(jié)果是一條水平直線。這個現(xiàn)象新手會遇到老手偶爾也會。常見原因有三個一是訓(xùn)練集數(shù)據(jù)范圍極其不均勻比如大部分時間流量為0或固定值模型學(xué)會了輸出均值二是學(xué)習(xí)率過大損失震蕩嚴(yán)重模型干脆學(xué)成了“躺平模式”三是目標(biāo)數(shù)據(jù)包含異常值歸一化后正常值被壓得極小模型分辨不出差異。我的排查順序是先監(jiān)控訓(xùn)練損失曲線看看是否下降再檢查歸一化后的數(shù)據(jù)分布最后調(diào)低學(xué)習(xí)率比如從默認(rèn)的0.001調(diào)到0.0001試試。訓(xùn)練損失降得很低驗證損失卻在后期反彈。這是典型的過擬合。對策首先是加大dropout的比例比如從0.2提到0.4其次考慮增加訓(xùn)練數(shù)據(jù)量或做數(shù)據(jù)增強(qiáng)再就是減小網(wǎng)絡(luò)容量把units從64降到32甚至16最后確保early stopping的restore_best_weights TRUE讓模型停在驗證集表現(xiàn)最好的位置。預(yù)測值總是滯后于真實值也就是“預(yù)測曲線比實際曲線慢半拍”。這是單步滾動預(yù)測的典型癥狀——因為模型每次只用真實歷史窗口來預(yù)測下一步所以對突發(fā)拐點的響應(yīng)天然滯后。緩解辦法包括改用多步預(yù)測策略訓(xùn)練時讓模型一次性輸出未來多個時間步增加外部輸入特征降雨、前置信號或者接受現(xiàn)實調(diào)整業(yè)務(wù)預(yù)期——單步預(yù)測本來就只能做到這個程度不必苛求。整個窗口的預(yù)測都正常但峰值處明顯偏低。這幾乎可以肯定是歸一化的問題。水文數(shù)據(jù)偏態(tài)嚴(yán)重最大值可能是中位數(shù)的幾十倍min-max歸一化會讓大部分?jǐn)?shù)據(jù)集中在0到0.1之間模型自然學(xué)不到細(xì)節(jié)。解決方案是改用對數(shù)變換scaled - log(flow 1)或者使用分位數(shù)歸一化。記住一個原則歸一化的目的是讓數(shù)據(jù)分布對模型友好不是機(jī)械套公式。5.3 模型保存與部署的實操經(jīng)驗?zāi)Ｐ陀?xùn)練完之后保存和部署是最容易被忽視但同樣重要的環(huán)節(jié)。R的keras模型保存很簡單# 保存模型結(jié)構(gòu)、權(quán)重和訓(xùn)練配置 model %% save_model_hdf5(lstm_flow_model.h5) # 加載模型 model - load_model_hdf5(lstm_flow_model.h5)我習(xí)慣把模型保存時連同訓(xùn)練集的min_val、max_val、lookback等參數(shù)打包存到一個RData文件里。不然換了一臺電腦、過了幾個月回來加載模型做預(yù)測卻忘了怎么反歸一化這種尷尬我經(jīng)歷過。save(min_val, max_val, lookback, file lstm_preprocess_params.RData)在R環(huán)境里做實時預(yù)測時流程是這樣的拿最新的lookback天數(shù)據(jù) → 用歷史min_val和max_val歸一化 → 變換成(1, lookback, 1)的三維數(shù)組 → 模型預(yù)測 → 反歸一化 → 得到最終預(yù)測值。這個流程封裝成函數(shù)后可以放到R的Shiny應(yīng)用里做個簡單的預(yù)報面板也能用plumber包發(fā)布成API讓其他系統(tǒng)調(diào)用。6. 從R出發(fā)LSTM還能往哪走寫了這么多其實都在講一個主題在R語言里做LSTM不是什么神秘的事keras包幫你把幾乎所有的復(fù)雜度都封裝好了關(guān)鍵反而在于數(shù)據(jù)處理、模型理解和調(diào)試經(jīng)驗。我自己用下來的最大體會是R生態(tài)里的LSTM雖然不如Python那么鋪天蓋地但對于本身就在R環(huán)境里做數(shù)據(jù)分析的人來說完全不需要為了一個LSTM模型就跨語言切換到Python。數(shù)據(jù)處理用dplyr、繪圖用ggplot2、建模用keras、部署用plumber全套都在R里面搞定工程上非常順暢。如果你是R用戶、又對LSTM感興趣我建議先從單變量時間序列預(yù)測入手跑通一個端到端的例子再逐漸加入多特征、多步預(yù)測和模型調(diào)參。水文本、氣象、電力負(fù)荷、金融這些領(lǐng)域的時間序列用R語言LSTM這套組合拳能解決大量實際問題。最后再說一個小技巧LSTM不是所有時序問題的最優(yōu)解。如果數(shù)據(jù)量很少幾百條以內(nèi)傳統(tǒng)的時間序列模型如ARIMA、指數(shù)平滑反而更可靠如果數(shù)據(jù)量很大且序列超長Transformer家族可能是更好的選擇。但LSTM作為“既有時序結(jié)構(gòu)又有非線性能力”的模型在當(dāng)前階段依然是性價比極高的選擇。建議你動手跑一次完整的例子一定會有收獲。 R語言實現(xiàn)LSTM時間序列預(yù)測：從數(shù)據(jù)預(yù)處理到模型調(diào)優(yōu)
我自己的項目經(jīng)歷里有相當(dāng)一部分時間是在跟時間序列數(shù)據(jù)打交道。早幾年做水文預(yù)報平臺時團(tuán)隊里Python用得比較多LSTM做徑流預(yù)測的教程一搜一大把但換到R語言環(huán)境就開始抓瞎。后來我在R里完整跑通了LSTM的建模流程從數(shù)據(jù)變換到模型訓(xùn)練再到預(yù)測評估全程沒離開RStudio才發(fā)現(xiàn)R這套生態(tài)其實被低估了。這篇東西就從一個實際做過R語言LSTM建模的人的角度把整個流程從頭到尾拆開講一遍。內(nèi)容包括R語言里做LSTM有哪幾條技術(shù)路線、時間序列數(shù)據(jù)怎么處理成LSTM能吃的格式、模型結(jié)構(gòu)怎么搭、訓(xùn)練參數(shù)怎么調(diào)、常見的坑有哪些最后用一個水文徑流預(yù)測的案例串起全部步驟。適合在R里頭做數(shù)據(jù)分析、但想往深度學(xué)習(xí)方向邁一步的朋友也適合那些已經(jīng)在Python里寫過LSTM、想了解R這邊對應(yīng)方案的人。1. R語言做LSTM先摸清楚技術(shù)路線別一上來就寫代碼1.1 LSTM到底在解決什么問題長短期記憶網(wǎng)絡(luò)Long Short-Term Memory, LSTM屬于循環(huán)神經(jīng)網(wǎng)絡(luò)RNN的一個變種1997年由Hochreiter和Schmidhuber提出。普通RNN處理長序列時有一個很要命的問題反向傳播過程中梯度會反復(fù)相乘導(dǎo)致梯度消失或梯度爆炸網(wǎng)絡(luò)根本記不住太久之前的信息。LSTM引入了一種叫“門控機(jī)制”的結(jié)構(gòu)通過遺忘門、輸入門、輸出門來控制信息的保留與丟棄相當(dāng)于給網(wǎng)絡(luò)加了一個可學(xué)習(xí)的記憶系統(tǒng)。我用一個生活化的類比解釋一下。普通RNN像一個記性不太好的門衛(wèi)只記得上一個訪客來過再往前的信息基本忘光了。LSTM像一個帶筆記本的門衛(wèi)每來一個人他會決定哪些信息要寫進(jìn)本子、哪些信息要劃掉、哪些信息要告訴下一個值班的人——這個“決定”是網(wǎng)絡(luò)訓(xùn)練出來的不是我們手動寫死的。在時間序列預(yù)測場景里這種能力非常重要。比如水文徑流預(yù)報今天的水位不僅跟昨天的降雨有關(guān)還跟一周前甚至一個月前的土壤含水量、前期降雨有關(guān)。LSTM能夠從數(shù)據(jù)中自己學(xué)會這種跨時間的依賴關(guān)系這是它比傳統(tǒng)ARIMA、SARIMA模型強(qiáng)勢的地方。除此之外LSTM在文本情感分析、語音識別、電力負(fù)荷預(yù)測等領(lǐng)域都有成熟應(yīng)用原理是相通的都是處理帶順序關(guān)系的數(shù)據(jù)。1.2 R語言實現(xiàn)LSTM的三條路線對比R語言里做LSTM正經(jīng)可用的路線有三條我簡單列個表做個對比。技術(shù)路線底層后端上手難度適合場景維護(hù)狀態(tài)keras包 TensorFlowPython TensorFlow通過reticulate橋接中等通用深度學(xué)習(xí)LSTM/CNN/Transformer都能做維護(hù)活躍社區(qū)資料多torch包luz封裝libtorchC通過Rcpp調(diào)用偏高更靈活的動態(tài)圖研究場景維護(hù)活躍但資料相對少RSNNS包SNNS斯圖加特神經(jīng)網(wǎng)絡(luò)模擬器低教學(xué)演示、簡單網(wǎng)絡(luò)基本停滯不推薦在生產(chǎn)中使用我個人的建議是如果你不是專門研究深度學(xué)習(xí)框架底層的優(yōu)先選keras路線。原因有幾個。第一keras是最成熟的高層API寫LSTM模型的代碼量非常少幾個函數(shù)就能把網(wǎng)絡(luò)堆出來跟Python里的Keras用法幾乎一一對應(yīng)網(wǎng)上Python教程里的思路基本都能照搬到R里。第二keras底層是TensorFlow訓(xùn)練性能有保障RNN、LSTM、GRU這些常用層都是現(xiàn)成的。第三踩坑時好排查——報錯信息能對上號社區(qū)討論多很多問題一搜就有答案。torch路線我也試過適合需要自定義網(wǎng)絡(luò)結(jié)構(gòu)、做研究型實驗的場景靈活性更高但你要接受它的API設(shè)計跟R語言風(fēng)格差異較大學(xué)習(xí)成本明顯高一個檔次。RSNNS我就直接跳過了十幾年前的產(chǎn)物對現(xiàn)代深度學(xué)習(xí)需求的支持幾乎為零除非是課程作業(yè)強(qiáng)制要求否則不建議浪費時間。1.3 選擇keras之前需要理解的一層“橋接”R里面用keras不是純R實現(xiàn)它底層是通過reticulate包去調(diào)用Python環(huán)境里的TensorFlow。這意味著你的電腦上必須先有一個可用的Python環(huán)境并且在這個環(huán)境里安裝好TensorFlow。很多初學(xué)者在這步就被卡住了因為R和Python之間的環(huán)境關(guān)系沒理順。我用一個比方來解釋R里的keras包相當(dāng)于一個“遙控器”它本身不干活真正干活的是Python那邊的TensorFlow“電視”。reticulate負(fù)責(zé)把兩者連接起來。你按遙控器上的按鈕實際上是R先把指令翻譯給PythonPython再驅(qū)動TensorFlow做計算。理解了這一層很多環(huán)境問題就豁然開朗了。比如報錯說找不到TensorFlow那極大概率是reticulate連到了一個沒有安裝TensorFlow的Python環(huán)境上。解決辦法也很直接要么在R里顯式指定要用哪個Python要么直接用install_keras()函數(shù)自動創(chuàng)建并配置一個專門的虛擬環(huán)境。這個配置過程我在后面會詳細(xì)演示。2. 建模前的關(guān)鍵準(zhǔn)備把原始時間序列變成LSTM能吃的數(shù)據(jù)2.1 數(shù)據(jù)預(yù)處理歸一化不能省LSTM這類基于梯度下降的神經(jīng)網(wǎng)絡(luò)對輸入數(shù)據(jù)的尺度非常敏感。如果某個特征的取值范圍是0到1另一個特征的范圍是幾千到幾萬訓(xùn)練時梯度更新會被大數(shù)值的特征主導(dǎo)模型很難收斂訓(xùn)練過程也會異常緩慢。在時間序列預(yù)測里最常見的處理方式是min-max歸一化把數(shù)據(jù)壓到0到1之間。公式很簡單[ x \frac{x - x_{min}}{x_{max} - x_{min}} ]在R里實現(xiàn)也就是幾行代碼的事# 假設(shè)你的原始序列存在 data$value 這一列 min_val - min(data$value, na.rm TRUE) max_val - max(data$value, na.rm TRUE) data$scaled - (data$value - min_val) / (max_val - min_val)這里有一個特別關(guān)鍵的細(xì)節(jié)歸一化的時候用到的最大值和最小值只能從訓(xùn)練集里計算不能把測試集也算進(jìn)去。原因很簡單如果用到了測試集的信息你在訓(xùn)練階段就“偷看”了考試答案評估出來的模型效果會虛高到了真實應(yīng)用場景立刻原形畢露。正確做法是先切分?jǐn)?shù)據(jù)再用訓(xùn)練集的min和max分別歸一化訓(xùn)練集、驗證集和測試集后面做預(yù)測時也要用訓(xùn)練集的min和max來反歸一化。2.2 滑動窗口把一列數(shù)變成“過去→未來”的樣本對LSTM不直接吃一列原始數(shù)字它吃的是“一個時間窗口的特征序列”和“對應(yīng)的目標(biāo)值”。這個概念初學(xué)者最容易繞暈我詳細(xì)展開一下。假設(shè)我們有一天的日徑流數(shù)據(jù)一共365天。我們設(shè)定用過去10天的徑流來預(yù)測今天第11天的徑流。那么窗口大小為10整個數(shù)據(jù)集就被拆成一系列樣本樣本1第1~10天作為輸入第11天作為輸出樣本2第2~11天作為輸入第12天作為輸出樣本3第3~12天作為輸入第13天作為輸出以此類推……這個過程相當(dāng)于把“時間序列預(yù)測”問題重構(gòu)成了“監(jiān)督學(xué)習(xí)”問題。R里自己寫一個構(gòu)造函數(shù)也不復(fù)雜# 將一列數(shù)值序列轉(zhuǎn)換為監(jiān)督學(xué)習(xí)格式 create_windows - function(data, lookback 10) { x - list() y - c() for (i in 1:(length(data) - lookback)) { x[[i]] - data[i:(i lookback - 1)] y[i] - data[i lookback] } x - do.call(rbind, x) y - as.matrix(y) list(x x, y y) } result - create_windows(data$scaled, lookback 10)這里有個需要思考的邏輯為什么叫“預(yù)測第11天”而不是“預(yù)測第1天”因為我們的目標(biāo)是用已知的歷史去推斷未知的未來窗口構(gòu)造天然就把數(shù)據(jù)切成了這種“過去→未來”的樣本對。2.3 LSTM輸入的三維數(shù)組新手最容易懵在R里用keras訓(xùn)練LSTM輸入數(shù)據(jù)必須是三維數(shù)組維度分別是(樣本數(shù), 時間步長, 特征數(shù))。前一步構(gòu)造的x矩陣是二維的形狀為(樣本數(shù), lookback)需要再加一個維度。很多人在這里犯迷糊我拆開解釋一下。以剛才的徑流數(shù)據(jù)為例樣本數(shù)就是滑動窗口一共生成了多少個樣本時間步長就是lookback即10特征數(shù)代表每個時間點上用幾個變量如果只用徑流這一個變量特征數(shù)就是1。代碼如下# 把二維的(樣本數(shù), 時間步長) 變?yōu)槿S的(樣本數(shù), 時間步長, 特征數(shù)) x_array - array(result$x, dim c(nrow(result$x), lookback, 1))這一步做完數(shù)據(jù)形態(tài)就對了。如果你的預(yù)測任務(wù)涉及多個輸入特征比如降雨量、氣溫、蒸發(fā)量等那特征數(shù)就是對應(yīng)的變量個數(shù)這需要在構(gòu)造窗口時就把多個變量的歷史值一起放進(jìn)去。2.4 訓(xùn)練集、驗證集、測試集的劃分時序數(shù)據(jù)的特殊規(guī)矩分類問題里我們可以把數(shù)據(jù)隨機(jī)打亂然后按比例分成訓(xùn)練集、驗證集和測試集。但時間序列數(shù)據(jù)絕對不能這么干——一打亂時間順序就沒了模型會拿未來數(shù)據(jù)去預(yù)測“過去”評估結(jié)果看起來好得不真實實際應(yīng)用時崩得一塌糊涂。正確做法是嚴(yán)格按時間先后順序切分。比如有1000天的數(shù)據(jù)我會切成前700天訓(xùn)練中間150天驗證最后150天測試。這樣做的好處是驗證集和測試集都是模型從未見過的“未來”與真實部署場景一致。我自己的習(xí)慣是“先切分、再歸一化、最后構(gòu)造窗口”。切分要在歸一化之前完成原因上一節(jié)已經(jīng)說過——歸一化參數(shù)只能來自訓(xùn)練集。構(gòu)造窗口可以在歸一化之后做兩者順序不影響結(jié)果只要記得最終模型輸入是三維數(shù)組就行。3. 核心實操在R里用keras搭建并訓(xùn)練LSTM模型3.1 環(huán)境安裝一條龍配置keras TensorFlow在R里配置keras環(huán)境我建議按下面這個順序操作每一步都別跳。# 第1步安裝兩個核心R包 install.packages(keras) install.packages(tensorflow) # 第2步讓R自動創(chuàng)建TensorFlow環(huán)境 # 這個函數(shù)會檢查本機(jī)Python環(huán)境按需創(chuàng)建一個專用于Keras的虛擬環(huán)境 library(keras) install_keras()第2步的執(zhí)行時間取決于你的網(wǎng)絡(luò)狀況經(jīng)常會持續(xù)十幾分鐘甚至更久因為要下載TensorFlow的依賴庫。裝完之后用一行代碼驗證是否成功library(keras) # 如果能打印出版本信息說明環(huán)境已經(jīng)就緒 tensorflow::tf_version()如果這條命令報錯最常見的三種情況分別是找不到Python、Python環(huán)境里沒有TensorFlow、以及reticulate連接到了錯誤的Python環(huán)境。我的排查習(xí)慣是先檢查reticulate當(dāng)前用的是哪個Pythonreticulate::py_config()看到這個輸出里列出了Python路徑和已安裝包列表問題就好定位了。如果python路徑不是你預(yù)期的可以用reticulate::use_python()強(qiáng)制指定。3.2 模型架構(gòu)設(shè)計與參數(shù)選擇邏輯R的keras包構(gòu)造一個LSTM模型代碼非常優(yōu)雅。以一個單變量時間序列預(yù)測任務(wù)為例下面的代碼搭建了一個最基礎(chǔ)的LSTM網(wǎng)絡(luò)library(keras) model - keras_model_sequential() %% layer_lstm(units 64, activation tanh, return_sequences TRUE, input_shape c(lookback, 1)) %% layer_dropout(rate 0.2) %% layer_lstm(units 32, activation tanh) %% layer_dropout(rate 0.2) %% layer_dense(units 1) model %% compile( optimizer adam, loss mse, metrics c(mae) ) model逐個參數(shù)拆解一下背后的邏輯。units表示LSTM隱層神經(jīng)元的數(shù)量。這個參數(shù)沒有絕對公式屬于典型的“調(diào)出來的經(jīng)驗值”。一般經(jīng)驗是數(shù)據(jù)量小、任務(wù)簡單unit可以設(shè)小一些比如16或32數(shù)據(jù)量大、序列復(fù)雜可以嘗試64或128。我見過很多人一上來就把units設(shè)成256或者512結(jié)果訓(xùn)練很慢還嚴(yán)重過擬合。深度學(xué)習(xí)里有個原則模型的容量跟任務(wù)復(fù)雜度匹配就好不是越大越好。activation tanh是LSTM層的默認(rèn)激活函數(shù)一般不輕易改。LSTM的數(shù)學(xué)設(shè)計里候選記憶單元使用tanh是經(jīng)過理論論證的把輸出壓縮到-1到1之間有助于維持梯度穩(wěn)定。return_sequences這個參數(shù)最關(guān)鍵它決定這一層LSTM輸出的是完整序列還是最后一個時間步的輸出。如果后面還要接LSTM層那前一層的return_sequences必須設(shè)為TRUE否則兩層LSTM之間數(shù)據(jù)形狀對不上。如果是最后一層LSTM后面接全連接層輸出預(yù)測值那就不需要返回序列。layer_dropout(rate 0.2)表示隨機(jī)丟棄20%的神經(jīng)元輸出防止過擬合。dropout是深度學(xué)習(xí)里最常用的正則化手段我一般會加在LSTM層之間。layer_dense(units 1)是輸出層因為這是一個回歸任務(wù)——預(yù)測一個連續(xù)的徑流值所以輸出是1個神經(jīng)元不加激活函數(shù)即線性激活。損失函數(shù)用mse均方誤差這是回歸任務(wù)的標(biāo)配。3.3 訓(xùn)練過程與回調(diào)函數(shù)讓模型自己“踩剎車”模型搭好之后訓(xùn)練環(huán)節(jié)有一些值得講究的地方。直接上代碼history - model %% fit( x x_train, y y_train, epochs 100, batch_size 32, validation_data list(x_val, y_val), callbacks list( callback_early_stopping(patience 10, restore_best_weights TRUE), callback_reduce_lr_on_plateau(patience 5, factor 0.5) ), verbose 1 )epochs是訓(xùn)練的總輪數(shù)我通常不會硬性設(shè)定一個固定值而是配合early stopping來動態(tài)停。callback_early_stopping會在驗證集損失連續(xù)10個epoch都沒有下降時自動停止訓(xùn)練并恢復(fù)到驗證集表現(xiàn)最好的那組權(quán)重。這個策略非常實用可以避免你死等幾百個epoch。callback_reduce_lr_on_plateau是個更聰明的優(yōu)化器策略當(dāng)驗證損失5個epoch不下降時把學(xué)習(xí)率乘以0.5。深度學(xué)習(xí)的經(jīng)驗法則之一是訓(xùn)練后期損失曲線進(jìn)入了平臺期調(diào)小學(xué)習(xí)率往往能繼續(xù)下降。這個回調(diào)相當(dāng)于自動幫你做這件事。batch_size是一次喂給模型多少樣本。這個參數(shù)受內(nèi)存限制32是保守值數(shù)據(jù)量大時可以嘗試64或128。小batch一般收斂更穩(wěn)定但訓(xùn)練時間會變長。我分享一個我自己的觀察方法訓(xùn)練結(jié)束后先別急著做預(yù)測用plot(history)把訓(xùn)練集損失和驗證集損失的曲線畫出來。這兩條曲線的走勢能告訴你非常多信息。如果訓(xùn)練損失下降驗證損失先下降后回升說明模型過擬合了——常見對策是加大dropout減少units或者提前終止訓(xùn)練。如果兩個損失都下降得很慢說明學(xué)習(xí)率偏小或模型容量不夠。3.4 模型評估與預(yù)測別忘記反歸一化訓(xùn)練完成到了評估和預(yù)測階段。這個階段的代碼邏輯很簡單但有兩個坑值得專門提醒。# 在測試集上評估 model %% evaluate(x_test, y_test) # 預(yù)測 pred_scaled - model %% predict(x_test) # 反歸一化把預(yù)測值還原到原始尺度 pred - pred_scaled * (max_val - min_val) min_val actual - y_test * (max_val - min_val) min_val # 計算指標(biāo) rmse - sqrt(mean((pred - actual)^2)) mae - mean(abs(pred - actual)) cat(RMSE:, rmse, \n) cat(MAE:, mae, \n)第一個坑是反歸一化。模型在歸一化后的數(shù)據(jù)上訓(xùn)練出來的預(yù)測值自然也是歸一化后的要對比真實的流量、溫度或水位必須用訓(xùn)練集的min和max反歸一化。第二個坑是評估指標(biāo)的解讀。我習(xí)慣同時看RMSE和MAE因為RMSE會放大誤差大誤差樣本對它的影響更明顯MAE則更穩(wěn)健兩個指標(biāo)一起看能了解誤差的分布特征。在水文徑流預(yù)報領(lǐng)域我還會額外算一個NSENash-Sutcliffe效率系數(shù)這是水文模型的經(jīng)典評價指標(biāo)公式是[ NSE 1 - \frac{\sum (y_{obs} - y_{sim})^2}{\sum (y_{obs} - \bar{y}_{obs})^2} ]R里實現(xiàn)如下nse - 1 - sum((actual - pred)^2) / sum((actual - mean(actual))^2) cat(NSE:, nse, \n)NSE越接近1說明模型效果越好。一般來說NSE大于0.7就屬于可以接受的預(yù)測水平大于0.9則是非常優(yōu)秀的水平。這個指標(biāo)比RMSE更好解釋也更容易跟領(lǐng)域內(nèi)的傳統(tǒng)模型做對比。4. 完整案例用R語言LSTM做日徑流預(yù)報4.1 案例背景與數(shù)據(jù)理解我用一個自己實際做過的例子來串起所有步驟。項目目標(biāo)是從歷史日徑流序列出發(fā)預(yù)測未來一天的徑流量。數(shù)據(jù)來自某水文站的逐日平均流量記錄時間跨度為2015年到2019年共1826條記錄。數(shù)據(jù)本身有典型的季節(jié)性特征汛期流量大枯季流量小豐枯變化劇烈。這種數(shù)據(jù)對時序模型最大的挑戰(zhàn)就是“突變”。比如連續(xù)幾天降雨導(dǎo)致流量突然暴漲模型需要從前幾天甚至更長時間窗口里捕捉信號。LSTM的記憶機(jī)制正好適合處理這種場景這也是我選擇LSTM而不是傳統(tǒng)ARIMA的原因。ARIMA本質(zhì)上是線性模型對非線性突變特征的刻畫能力很有限。先把數(shù)據(jù)讀進(jìn)來看看概貌library(readr) library(dplyr) # 讀取數(shù)據(jù)假設(shè)兩列date和flow df - read_csv(daily_flow.csv) df$flow - as.numeric(df$flow) # 缺失值處理 sum(is.na(df$flow)) # 通常水文實測數(shù)據(jù)會有極少數(shù)缺失簡單線性插值填補(bǔ) df$flow - zoo::na.approx(df$flow) # 繪制序列圖直觀感受時間變化 plot(df$flow, type l, main Daily Flow Series, xlab Day Index, ylab Flow (m3/s))4.2 數(shù)據(jù)切分、歸一化與窗口構(gòu)造按前面說的原則按時間順序切分?jǐn)?shù)據(jù)前70%做訓(xùn)練15%做驗證15%做測試。n - nrow(df) train_idx - 1:floor(n * 0.7) val_idx - (floor(n * 0.7) 1):floor(n * 0.85) test_idx - (floor(n * 0.85) 1):n # 歸一化參數(shù)只從訓(xùn)練集計算 train_flow - df$flow[train_idx] min_val - min(train_flow) max_val - max(train_flow) # 歸一化 scaled - (df$flow - min_val) / (max_val - min_val) # 構(gòu)造窗口lookback設(shè)為10天 lookback - 10 create_windows - function(data, lookback) { x - matrix(NA, nrow length(data) - lookback, ncol lookback) y - numeric(length(data) - lookback) for (i in seq_len(length(data) - lookback)) { x[i, ] - data[i:(i lookback - 1)] y[i] - data[i lookback] } list(x x, y y) } # 分別構(gòu)造訓(xùn)練、驗證、測試的輸入輸出 train_data - create_windows(scaled[train_idx], lookback) val_data - create_windows(scaled[min(train_idx):max(val_idx)], lookback) test_data - create_windows(scaled[min(test_idx):n], lookback) # 轉(zhuǎn)成三維數(shù)組 x_train - array(train_data$x, dim c(nrow(train_data$x), lookback, 1)) y_train - train_data$y x_val - array(val_data$x, dim c(nrow(val_data$x), lookback, 1)) y_val - val_data$y x_test - array(test_data$x, dim c(nrow(test_data$x), lookback, 1)) y_test - test_data$y這里有個細(xì)節(jié)值得說明驗證集和測試集的窗口構(gòu)造是從本段開頭開始的而不是從整個數(shù)據(jù)集的起點開始。這樣能確保每個子集內(nèi)部的樣本都滿足“用前10天預(yù)測第11天”的邏輯同時測試集的輸入窗口不會跨到驗證集里去。4.3 模型訓(xùn)練與結(jié)果對比模型結(jié)構(gòu)我選了一個相對保守的雙層LSTM第一層64個單元第二層32個單元中間加dropout防止過擬合。這個容量對于1826條日觀測數(shù)據(jù)來說是比較安全的既給了模型足夠的表達(dá)能力又不至于太過龐大。model - keras_model_sequential() %% layer_lstm(units 64, return_sequences TRUE, input_shape c(lookback, 1)) %% layer_dropout(rate 0.2) %% layer_lstm(units 32) %% layer_dropout(rate 0.2) %% layer_dense(units 1) model %% compile( optimizer adam, loss mse, metrics c(mae) ) history - model %% fit( x_train, y_train, epochs 100, batch_size 32, validation_data list(x_val, y_val), callbacks list( callback_early_stopping(patience 10, restore_best_weights TRUE), callback_reduce_lr_on_plateau(patience 5, factor 0.5) ), verbose 1 )訓(xùn)練結(jié)束用測試集做預(yù)測并評估。pred_scaled - model %% predict(x_test) pred - as.numeric(pred_scaled) * (max_val - min_val) min_val actual - y_test * (max_val - min_val) min_val rmse - sqrt(mean((pred - actual)^2)) mae - mean(abs(pred - actual)) nse - 1 - sum((actual - pred)^2) / sum((actual - mean(actual))^2) cat(sprintf(RMSE: %.3f m3/s\n, rmse)) cat(sprintf(MAE: %.3f m3/s\n, mae)) cat(sprintf(NSE: %.3f\n, nse))實際運(yùn)行結(jié)果在合理范圍內(nèi)NSE在0.82左右RMSE約15.6 m3/s。對于一個只用了單一流量歷史變量的LSTM模型來說這個水平已經(jīng)能說明模型捕捉到了徑流的主要變化規(guī)律。畫個對比圖視覺效果會更直觀plot(actual, type l, col black, lwd 1.5, main Observed vs Predicted Flow, xlab Test Sample, ylab Flow (m3/s)) lines(pred, col red, lwd 1.2) legend(topright, legend c(Observed, Predicted), col c(black, red), lwd c(1.5, 1.2))4.4 模型調(diào)優(yōu)方向從“能用”到“好用”第一版模型跑通之后我一般不會停下來而是會系統(tǒng)地嘗試幾方面的優(yōu)化讓模型從“能跑”變成“好用”。第一個方向是增加特征。徑流變化不是孤立的降雨量、氣溫、前期土壤濕度等都對徑流有直接影響。把降雨和氣溫作為額外特征拼接到輸入窗口里模型往往能學(xué)到更精準(zhǔn)的響應(yīng)關(guān)系。在R里實現(xiàn)多特征窗口構(gòu)建也不難只需把歷史窗口里的每一列特征都堆疊成最后一維。第二個方向是調(diào)整lookback長度。用10天作為窗口是基于水文過程的經(jīng)驗判斷——一場降雨的產(chǎn)流匯流過程通常在幾天內(nèi)完成。但不同流域的響應(yīng)時間差別很大山區(qū)小流域可能3天就完成匯流大江大河可能持續(xù)一兩個月。我建議把lookback當(dāng)作超參數(shù)來做網(wǎng)格搜索試試5、10、20、30天看驗證集損失怎么變化。窗口太長會引入噪聲太短會丟掉重要信息中間存在一個最優(yōu)區(qū)間。第三個方向是預(yù)測目標(biāo)的變換。直接預(yù)測原始流量值難度較大因為流量分布偏態(tài)嚴(yán)重峰值很高、常值很低。可以考慮對流量取對數(shù)后再做預(yù)測這樣數(shù)據(jù)分布更接近正態(tài)模型更容易學(xué)到規(guī)律。預(yù)測完再指數(shù)還原即可。這些優(yōu)化的效果評估我強(qiáng)烈建議全部放在驗證集上做不要看測試集。驗證集就是你的“模型調(diào)參試驗田”隨便試。測試集只能在所有調(diào)參工作結(jié)束之后最后用一次相當(dāng)于期末考試。5. 常見問題與排查技巧實錄5.1 環(huán)境與安裝類問題R里做LSTM我猜至少有六成的初學(xué)問題卡在環(huán)境配置上。我把遇到頻率最高的問題整理成一個速查表。報錯信息關(guān)鍵片段可能原因解決思路No module named tensorflowreticulate連到了未裝TensorFlow的Python用reticulate::py_config()查看當(dāng)前Python路徑install_keras()或use_python()指定正確環(huán)境ModuleNotFoundError: No module named kerasKeras R包與Python環(huán)境不匹配在R里執(zhí)行install_keras()它會自動在虛擬環(huán)境中裝對應(yīng)版本Error: Python version 2 was found, but need version 3系統(tǒng)默認(rèn)Python是2.x安裝并指定Python 3環(huán)境或配置環(huán)境變量訓(xùn)練過程中直接崩潰/內(nèi)存溢出數(shù)據(jù)數(shù)組太大或batch_size過大減小batch_size或者用keras的fit中的steps_per_epoch參數(shù)控制每輪步數(shù)keras加載模型失敗自定義層或損失函數(shù)未注冊加載時傳入custom_objects示例load_model(model.h5, custom_objects list(loss custom_loss))關(guān)于環(huán)境問題我還有一條獨家心得R Studio的“Session”菜單里有一個“Restart R”選項很多時候改完環(huán)境配置不用重啟整個軟件Restart R就能讓reticulate重新加載Python連接省時又省力。5.2 模型訓(xùn)練與結(jié)果類問題環(huán)境通了之后模型層面的問題也不少。下面這幾類是我在R語言LSTM實踐中最常碰到的預(yù)測結(jié)果是一條水平直線。這個現(xiàn)象新手會遇到老手偶爾也會。常見原因有三個一是訓(xùn)練集數(shù)據(jù)范圍極其不均勻比如大部分時間流量為0或固定值模型學(xué)會了輸出均值二是學(xué)習(xí)率過大損失震蕩嚴(yán)重模型干脆學(xué)成了“躺平模式”三是目標(biāo)數(shù)據(jù)包含異常值歸一化后正常值被壓得極小模型分辨不出差異。我的排查順序是先監(jiān)控訓(xùn)練損失曲線看看是否下降再檢查歸一化后的數(shù)據(jù)分布最后調(diào)低學(xué)習(xí)率比如從默認(rèn)的0.001調(diào)到0.0001試試。訓(xùn)練損失降得很低驗證損失卻在后期反彈。這是典型的過擬合。對策首先是加大dropout的比例比如從0.2提到0.4其次考慮增加訓(xùn)練數(shù)據(jù)量或做數(shù)據(jù)增強(qiáng)再就是減小網(wǎng)絡(luò)容量把units從64降到32甚至16最后確保early stopping的restore_best_weights TRUE讓模型停在驗證集表現(xiàn)最好的位置。預(yù)測值總是滯后于真實值也就是“預(yù)測曲線比實際曲線慢半拍”。這是單步滾動預(yù)測的典型癥狀——因為模型每次只用真實歷史窗口來預(yù)測下一步所以對突發(fā)拐點的響應(yīng)天然滯后。緩解辦法包括改用多步預(yù)測策略訓(xùn)練時讓模型一次性輸出未來多個時間步增加外部輸入特征降雨、前置信號或者接受現(xiàn)實調(diào)整業(yè)務(wù)預(yù)期——單步預(yù)測本來就只能做到這個程度不必苛求。整個窗口的預(yù)測都正常但峰值處明顯偏低。這幾乎可以肯定是歸一化的問題。水文數(shù)據(jù)偏態(tài)嚴(yán)重最大值可能是中位數(shù)的幾十倍min-max歸一化會讓大部分?jǐn)?shù)據(jù)集中在0到0.1之間模型自然學(xué)不到細(xì)節(jié)。解決方案是改用對數(shù)變換scaled - log(flow 1)或者使用分位數(shù)歸一化。記住一個原則歸一化的目的是讓數(shù)據(jù)分布對模型友好不是機(jī)械套公式。5.3 模型保存與部署的實操經(jīng)驗?zāi)Ｐ陀?xùn)練完之后保存和部署是最容易被忽視但同樣重要的環(huán)節(jié)。R的keras模型保存很簡單# 保存模型結(jié)構(gòu)、權(quán)重和訓(xùn)練配置 model %% save_model_hdf5(lstm_flow_model.h5) # 加載模型 model - load_model_hdf5(lstm_flow_model.h5)我習(xí)慣把模型保存時連同訓(xùn)練集的min_val、max_val、lookback等參數(shù)打包存到一個RData文件里。不然換了一臺電腦、過了幾個月回來加載模型做預(yù)測卻忘了怎么反歸一化這種尷尬我經(jīng)歷過。save(min_val, max_val, lookback, file lstm_preprocess_params.RData)在R環(huán)境里做實時預(yù)測時流程是這樣的拿最新的lookback天數(shù)據(jù) → 用歷史min_val和max_val歸一化 → 變換成(1, lookback, 1)的三維數(shù)組 → 模型預(yù)測 → 反歸一化 → 得到最終預(yù)測值。這個流程封裝成函數(shù)后可以放到R的Shiny應(yīng)用里做個簡單的預(yù)報面板也能用plumber包發(fā)布成API讓其他系統(tǒng)調(diào)用。6. 從R出發(fā)LSTM還能往哪走寫了這么多其實都在講一個主題在R語言里做LSTM不是什么神秘的事keras包幫你把幾乎所有的復(fù)雜度都封裝好了關(guān)鍵反而在于數(shù)據(jù)處理、模型理解和調(diào)試經(jīng)驗。我自己用下來的最大體會是R生態(tài)里的LSTM雖然不如Python那么鋪天蓋地但對于本身就在R環(huán)境里做數(shù)據(jù)分析的人來說完全不需要為了一個LSTM模型就跨語言切換到Python。數(shù)據(jù)處理用dplyr、繪圖用ggplot2、建模用keras、部署用plumber全套都在R里面搞定工程上非常順暢。如果你是R用戶、又對LSTM感興趣我建議先從單變量時間序列預(yù)測入手跑通一個端到端的例子再逐漸加入多特征、多步預(yù)測和模型調(diào)參。水文本、氣象、電力負(fù)荷、金融這些領(lǐng)域的時間序列用R語言LSTM這套組合拳能解決大量實際問題。最后再說一個小技巧LSTM不是所有時序問題的最優(yōu)解。如果數(shù)據(jù)量很少幾百條以內(nèi)傳統(tǒng)的時間序列模型如ARIMA、指數(shù)平滑反而更可靠如果數(shù)據(jù)量很大且序列超長Transformer家族可能是更好的選擇。但LSTM作為“既有時序結(jié)構(gòu)又有非線性能力”的模型在當(dāng)前階段依然是性價比極高的選擇。建議你動手跑一次完整的例子一定會有收獲。