簡介這份資源面向具備一定R語言基礎、希望深入掌握時間序列建模的數(shù)據(jù)分析與預測從業(yè)者聚焦加法與乘法過程兩類核心思路并延伸至廣義可加模型GAM的非線性建模場景。壓縮包內共1個文件為R腳本整體約1KB可直接在R環(huán)境中加載運行便于對照代碼理解建模流程。內容圍繞時間序列的基本構成展開涉及趨勢、季節(jié)性與隨機成分的線性相加與相互影響兩種假設并串聯(lián)ARIMA、SARIMA、STL分解以及基于mgcv包的gam()擬合等實現(xiàn)路徑同時涵蓋數(shù)據(jù)加載、模型選擇、參數(shù)調優(yōu)與殘差診斷等環(huán)節(jié)。目前已有618人學習下載適合作為經(jīng)濟、金融、氣象等領域預測任務的實操參考幫助讀者把抽象模型轉化為可復用的代碼框架并借助GAM樣條項捕捉復雜非線性趨勢。1. 時間序列的加法與乘法為什么你的預測總在趨勢拐點翻車很多人做時間序列預測時拿到數(shù)據(jù)第一件事就是往 ARIMA 或 Prophet 里一塞跑出來一條平滑曲線看著擬合得不錯一到真實業(yè)務里就發(fā)現(xiàn)趨勢一變預測全廢。問題往往不在模型本身而在于你默認了一個假設——這個序列是「加法結構」還是「乘法結構」。這兩個詞聽起來像統(tǒng)計學課本里的老古董但在 R 語言里用 GAM廣義可加模型做時間序列分解和預測時它直接決定了你的模型是「能跟著趨勢走」還是「一到旺季就低估」。加法模型假設各成分獨立疊加觀測值 趨勢 季節(jié) 殘差。乘法模型假設成分之間是比例關系觀測值 趨勢 × 季節(jié) × 殘差。現(xiàn)實業(yè)務里銷售額隨趨勢增長時季節(jié)性波動幅度也在放大這就是典型的乘法結構。如果你硬用加法模型去擬合殘差會呈現(xiàn)明顯的異方差預測區(qū)間在高峰期窄得離譜低谷期又寬得沒用。R 語言里的 GAM 通過mgcv包提供了靈活的平滑項可以讓你在同一個框架下處理這兩種結構甚至混合結構。這篇文章面向的是已經(jīng)會用 R 做基本數(shù)據(jù)分析、但想在時間序列上把 GAM 用對的從業(yè)者。我會把加法與乘法的選擇邏輯、GAM 的平滑項設置、以及實際跑通一套分解加預測的流程講清楚順帶把那些讓我翻過車的參數(shù)坑標出來。2. 加法與乘法結構的判定先看殘差再談模型2.1 從業(yè)務場景反推結構類型拿到一條時間序列不要急著畫圖。先問一個業(yè)務問題當趨勢上升 10% 時季節(jié)波動的絕對幅度是保持不變還是也跟著漲 10%如果是前者加法結構更合理如果是后者乘法結構更合理。這個判斷不需要任何統(tǒng)計檢驗靠業(yè)務常識就能定調。舉個例子某電商平臺的日訂單量。過去兩年整體趨勢從每天 1000 單漲到 3000 單同時「雙十一」當天的峰值從 5000 單漲到 15000 單。峰值與均值的比例大致穩(wěn)定在 5 倍左右這就是乘法結構的典型信號。反過來如果是一個成熟產(chǎn)品的日活用戶趨勢基本平緩周末比工作日固定少 2000 人那加法結構就夠了。在 R 里你可以用forecast包的decompose()函數(shù)分別跑加法和乘法分解然后對比殘差圖。乘法分解的殘差如果比加法分解更接近白噪聲且殘差絕對值不隨趨勢增大而膨脹那就選乘法。library(forecast) # 假設 ts_data 是一個 ts 對象頻率為 7周季節(jié) # 加法分解 decomp_add - decompose(ts_data, type additive) # 乘法分解 decomp_mult - decompose(ts_data, type multiplicative) # 對比殘差的標準差隨趨勢的變化 # 如果乘法分解的殘差更穩(wěn)定選乘法 sd_add - sd(decomp_add$random, na.rm TRUE) sd_mult - sd(decomp_mult$random, na.rm TRUE) # 更直觀的方法看殘差與趨勢的相關系數(shù) trend_add - decomp_add$trend resid_add - decomp_add$random cor_add - cor(trend_add, abs(resid_add), use complete.obs) trend_mult - decomp_mult$trend resid_mult - decomp_mult$random cor_mult - cor(trend_mult, abs(resid_mult), use complete.obs) cat(加法殘差與趨勢相關系數(shù):, cor_add, \n) cat(乘法殘差與趨勢相關系數(shù):, cor_mult, \n)這段代碼的邏輯是如果殘差的絕對值和趨勢有正相關說明波動幅度隨趨勢增長加法模型就不合適。參數(shù)上type指定分解類型frequency在ts()里設定。注意decompose()只適合季節(jié)性固定的序列如果季節(jié)模式在變得用stl()或 GAM。2.2 用 GAM 的平滑項同時捕捉趨勢與季節(jié)GAM 的優(yōu)勢在于它不預設趨勢是線性的也不預設季節(jié)是固定的。mgcv包里的gam()函數(shù)可以用s()指定平滑項用te()指定交互項。對于時間序列我一般這樣建library(mgcv) # 構造時間索引和季節(jié)索引 n - length(ts_data) time_idx - 1:n season_idx - cycle(ts_data) # 1 到 frequency # 加法結構 GAM gam_add - gam(ts_data ~ s(time_idx, k 20) s(season_idx, k 7, bs cc), family gaussian()) # 乘法結構先取對數(shù)再跑加法 GAM gam_mult - gam(log(ts_data) ~ s(time_idx, k 20) s(season_idx, k 7, bs cc), family gaussian())這里的關鍵參數(shù)是k它控制平滑項的最大自由度。k太小會欠擬合趨勢被抹平k太大會過擬合把噪聲當信號。我一般從 10 到 30 之間試用gam.check()看殘差和k的顯著性。bs cc是循環(huán)平滑專門用于季節(jié)索引這種周期性變量不加這個12 月和 1 月之間會出現(xiàn)斷點。乘法結構取對數(shù)后跑加法 GAM等價于在原始尺度上做乘法分解。預測時記得exp()回來但要注意偏差修正——直接exp()會低估均值因為對數(shù)正態(tài)分布的均值是exp(mu sigma^2/2)。這個坑我后面會細說。2.3 模型選擇的量化依據(jù)AIC 與殘差診斷加法 GAM 和乘法 GAM 跑完后不能只看圖。AIC()可以比較但注意乘法模型是在對數(shù)尺度上算的AIC 不能直接和加法模型比。正確做法是把兩個模型的預測值都變換回原始尺度算 RMSE 或 MAE用交叉驗證比較。# 樣本外交叉驗證留出最后 30 個點 train_n - n - 30 train_data - ts_data[1:train_n] test_data - ts_data[(train_n 1):n] # 重新擬合加法模型 gam_add_train - gam(train_data ~ s(1:train_n, k 20) s(cycle(train_data), k 7, bs cc)) # 重新擬合乘法模型 gam_mult_train - gam(log(train_data) ~ s(1:train_n, k 20) s(cycle(train_data), k 7, bs cc)) # 預測 pred_add - predict(gam_add_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(test_data) )) pred_mult_log - predict(gam_mult_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(test_data) )) pred_mult - exp(pred_mult_log) # 計算 RMSE rmse_add - sqrt(mean((test_data - pred_add)^2)) rmse_mult - sqrt(mean((test_data - pred_mult)^2)) cat(加法 RMSE:, rmse_add, \n) cat(乘法 RMSE:, rmse_mult, \n)這段代碼的核心是「用樣本外誤差說話」。參數(shù)上k在訓練集上重新選不要用全量數(shù)據(jù)調好的k直接套。cycle()提取季節(jié)位置predict()的newdata必須包含和訓練時同名的變量。如果乘法模型的 RMSE 明顯更低且殘差沒有異方差那就選乘法。3. 在 R 里跑通 GAM 時間序列從數(shù)據(jù)到預測的完整鏈路3.1 數(shù)據(jù)準備與 ts 對象構造R 里做時間序列第一步是把數(shù)據(jù)轉成ts對象。很多人從 CSV 讀進來直接跑結果cycle()返回 NULL季節(jié)平滑項直接報錯。正確做法是明確頻率和起始時間。# 假設 raw_data 是數(shù)據(jù)框date 列是日期value 列是觀測值 raw_data$date - as.Date(raw_data$date) # 按日期排序 raw_data - raw_data[order(raw_data$date), ] # 構造 ts 對象頻率為 7周數(shù)據(jù) ts_data - ts(raw_data$value, frequency 7, start c(year(min(raw_data$date)), as.numeric(format(min(raw_data$date), %j)))) # 檢查 cat(頻率:, frequency(ts_data), \n) cat(周期數(shù):, length(ts_data) / frequency(ts_data), \n) head(cycle(ts_data))frequency 7表示每周 7 個觀測start參數(shù)指定起始年份和年內第幾天。如果頻率設錯比如日數(shù)據(jù)設成 365cycle()會返回 1 到 365季節(jié)平滑項k就得設得很大計算量爆炸且容易過擬合。常見做法是日數(shù)據(jù)用frequency 7捕捉周內模式或者frequency 365.25捕捉年內模式但后者需要k至少 50 以上我一般先用周頻率跑通再考慮年頻率。3.2 GAM 平滑項的參數(shù)設置與 gam.check 解讀gam()函數(shù)里最關鍵的三個參數(shù)k、bs、family。k是基函數(shù)維度決定平滑曲線的最大彎曲次數(shù)。bs是基函數(shù)類型tp是薄板回歸樣條默認cc是循環(huán)三次樣條用于周期變量cr是三次回歸樣條計算更快。family指定分布族高斯用于連續(xù)值泊松用于計數(shù)負二項用于過離散計數(shù)。# 完整模型擬合 gam_fit - gam(ts_data ~ s(time_idx, k 25, bs tp) s(season_idx, k 7, bs cc), family gaussian(), method REML) # 用 REML 估計平滑參數(shù) # 檢查 gam.check(gam_fit)gam.check()輸出四張圖殘差 vs 擬合值、QQ 圖、直方圖、響應 vs 擬合值。重點看第一張如果殘差呈現(xiàn)漏斗形隨擬合值增大而擴散說明方差非恒定加法高斯模型不合適要么換乘法取對數(shù)要么換family Gamma。QQ 圖如果尾部偏離說明殘差非正態(tài)預測區(qū)間會不準。method REML是我強烈建議的默認的 GCV 在樣本量小的時候容易過擬合REML 更穩(wěn)健。這個參數(shù)不設gam()會用 GCV殘差診斷經(jīng)常顯示k不夠其實換了 REML 就好了。3.3 預測與置信區(qū)間對數(shù)變換后的偏差修正乘法模型預測時predict()返回的是對數(shù)尺度上的值exp()回去之后得到的是中位數(shù)不是均值。如果業(yè)務要的是均值預測必須加偏差修正。# 乘法模型預測對數(shù)尺度 pred_log - predict(gam_mult, newdata data.frame( time_idx (n 1):(n 30), season_idx rep(1:7, length.out 30) ), se.fit TRUE) # 偏差修正exp(mu sigma^2/2) sigma2 - summary(gam_mult)$scale pred_mean - exp(pred_log$fit sigma2 / 2) # 置信區(qū)間近似 lower - exp(pred_log$fit - 1.96 * pred_log$se.fit) upper - exp(pred_log$fit 1.96 * pred_log$se.fit) # 注意這個區(qū)間是條件均值區(qū)間不是預測區(qū)間 # 預測區(qū)間還要加殘差方差 pred_lower - exp(pred_log$fit - 1.96 * sqrt(pred_log$se.fit^2 sigma2)) pred_upper - exp(pred_log$fit 1.96 * sqrt(pred_log$se.fit^2 sigma2))summary(gam_mult)$scale提取的是殘差方差估計。偏差修正項sigma2 / 2看起來小但當sigma2大的時候不修正會低估均值 5% 到 10%在庫存補貨場景里就是真金白銀的差距。置信區(qū)間和預測區(qū)間是兩回事前者是均值的區(qū)間后者是單次觀測的區(qū)間。業(yè)務上要「明天可能賣多少」用預測區(qū)間要「明天平均賣多少」用置信區(qū)間。4. 避坑與排查那些讓我重跑模型的參數(shù)陷阱4.1 現(xiàn)象gam.check 顯示 k 不夠調大后過擬合原因k的默認值是 10對于趨勢變化劇烈的序列10 個基函數(shù)不夠彎曲。但直接調到 50模型會把噪聲也擬合進去樣本外誤差反而上升。解決先用gam.check()看殘差是否還有模式。如果殘差已經(jīng)像白噪聲k就夠了不用管 p 值。如果殘差還有趨勢每次加 5直到殘差無模式。同時用method REML它比 GCV 更不容易過擬合。我一般從k 20開始趨勢項和季節(jié)項分開調。4.2 現(xiàn)象乘法模型預測值在低谷期為負原因對數(shù)變換后預測再exp()回去理論上不會負。但如果用了family gaussian()且沒有取對數(shù)直接跑乘法結構比如把季節(jié)項寫成乘積形式線性預測子可能為負。解決乘法結構必須在對數(shù)尺度上建模。log(ts_data)后跑加法 GAMexp()回來。如果原始數(shù)據(jù)有零或負值先加一個常數(shù)再取對數(shù)或者改用family Gamma(link log)后者直接在對數(shù)鏈接函數(shù)上建模不需要手動變換。4.3 現(xiàn)象季節(jié)平滑項 bs cc 報錯「A term has fewer unique covariate combinations than specified maximum degrees of freedom」原因k設得比季節(jié)周期的唯一值數(shù)量還大。比如周頻率frequency 7season_idx只有 7 個唯一值k設成 10 就報錯。解決k必須小于等于唯一值數(shù)量。周頻率用k 7月頻率用k 12日頻率用k 7周內模式或k 365年內模式但計算量大。如果確實需要更大的k改用bs tp并手動構造傅里葉項。4.4 現(xiàn)象預測區(qū)間在趨勢上升段明顯偏窄原因GAM 的predict(se.fit TRUE)只給了平滑項的不確定性沒有包含殘差方差。而且如果模型是加法結構但數(shù)據(jù)實際是乘法結構殘差方差會隨趨勢增大區(qū)間估計系統(tǒng)性偏窄。解決預測區(qū)間用sqrt(se.fit^2 sigma2)不要只用se.fit。同時檢查殘差是否異方差如果是換乘法結構或family Gamma。我習慣在預測后畫殘差 vs 擬合值圖確認沒有漏斗形才發(fā)結果。4.5 現(xiàn)象時間索引 time_idx 從 1 到 n預測時新數(shù)據(jù)的時間索引對不上原因訓練時time_idx是 1 到 n預測時新數(shù)據(jù)的time_idx必須從 n1 開始。如果直接用1:30模型會以為你在預測歷史區(qū)間。解決預測時構造newdata的time_idx用(n1):(nh)season_idx用cycle()或手動指定。更穩(wěn)妥的做法是把時間索引存成數(shù)據(jù)框的一列訓練和預測用同一套構造邏輯。5. 進階技巧用 te() 捕捉趨勢與季節(jié)的交互效應5.1 什么時候需要交互項加法 GAM 假設趨勢和季節(jié)是獨立作用的趨勢上升不影響季節(jié)波動的形狀。但現(xiàn)實中很多序列的季節(jié)模式會隨趨勢變化。比如一個產(chǎn)品剛上市時只有周末賣得好隨著市場成熟工作日銷量也上來了季節(jié)波動的幅度和形狀都變了。這時候s(time_idx) s(season_idx)就不夠了需要用te(time_idx, season_idx)張量積平滑來捕捉交互。# 張量積交互模型 gam_te - gam(ts_data ~ te(time_idx, season_idx, k c(15, 7), bs c(tp, cc)), family gaussian(), method REML) # 對比無交互模型 gam_no_te - gam(ts_data ~ s(time_idx, k 15) s(season_idx, k 7, bs cc), family gaussian(), method REML) # 用 AIC 比較同一尺度下 cat(交互模型 AIC:, AIC(gam_te), \n) cat(無交互 AIC:, AIC(gam_no_te), \n)te()的參數(shù)k是一個向量分別對應兩個變量的基函數(shù)維度。bs也是向量對應各自的基函數(shù)類型。交互模型的自由度是k[1] * k[2]的量級計算量比加法模型大很多樣本量少于 200 時慎用容易過擬合。我一般先跑加法模型如果殘差在特定時間段比如旺季還有系統(tǒng)偏差再上交互項。5.2 用 vis.gam 可視化交互效應mgcv自帶的vis.gam()可以畫三維透視圖或等高線圖直觀看到季節(jié)模式如何隨趨勢變化。# 三維透視圖 vis.gam(gam_te, view c(time_idx, season_idx), theta 30, phi 20, color heat, plot.type persp, ticktype detailed) # 等高線圖 vis.gam(gam_te, view c(time_idx, season_idx), plot.type contour, color heat)view指定兩個變量theta和phi控制視角plot.type選persp或contour。等高線圖里如果等高線是平行的說明沒有交互如果彎曲或交叉說明季節(jié)模式隨趨勢變了。這個圖我每次做完交互模型都會看比 AIC 更直觀。5.3 一個我常犯的錯誤交互項加了但沒做預測對比剛用te()的時候我看到 AIC 降了就以為交互模型更好直接拿去預測結果樣本外 RMSE 反而更高。后來養(yǎng)成習慣不管 AIC 多低一定做樣本外交叉驗證。交互模型的方差更大樣本量不夠時AIC 的懲罰項不足以抵消過擬合。# 樣本外對比 train_n - n - 30 gam_te_train - gam(ts_data[1:train_n] ~ te(1:train_n, cycle(ts_data[1:train_n]), k c(15, 7), bs c(tp, cc)), family gaussian(), method REML) pred_te - predict(gam_te_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(ts_data[(train_n 1):n]) )) rmse_te - sqrt(mean((ts_data[(train_n 1):n] - pred_te)^2)) cat(交互模型樣本外 RMSE:, rmse_te, \n)如果rmse_te比加法模型的樣本外 RMSE 還大說明交互項在擬合噪聲果斷回退。這個習慣讓我少了很多「AIC 好看但業(yè)務翻車」的情況。5.4 最后的習慣先畫圖再跑模型我現(xiàn)在拿到任何時間序列第一件事是plot(ts_data)第二件事是plot(decompose(ts_data))或plot(stl(ts_data, s.window periodic))。圖上看一眼趨勢和季節(jié)的形態(tài)比任何統(tǒng)計檢驗都快。加法還是乘法很多時候圖上一眼就能定如果季節(jié)波動的幅度隨趨勢明顯放大乘法如果幅度穩(wěn)定加法。GAM 只是把這個判斷量化并給出預測區(qū)間。希望幫到你。本文還有配套的精品資源點擊獲取 R語言GAM時間序列預測：加法與乘法結構選擇及避坑指南
簡介這份資源面向具備一定R語言基礎、希望深入掌握時間序列建模的數(shù)據(jù)分析與預測從業(yè)者聚焦加法與乘法過程兩類核心思路并延伸至廣義可加模型GAM的非線性建模場景。壓縮包內共1個文件為R腳本整體約1KB可直接在R環(huán)境中加載運行便于對照代碼理解建模流程。內容圍繞時間序列的基本構成展開涉及趨勢、季節(jié)性與隨機成分的線性相加與相互影響兩種假設并串聯(lián)ARIMA、SARIMA、STL分解以及基于mgcv包的gam()擬合等實現(xiàn)路徑同時涵蓋數(shù)據(jù)加載、模型選擇、參數(shù)調優(yōu)與殘差診斷等環(huán)節(jié)。目前已有618人學習下載適合作為經(jīng)濟、金融、氣象等領域預測任務的實操參考幫助讀者把抽象模型轉化為可復用的代碼框架并借助GAM樣條項捕捉復雜非線性趨勢。1. 時間序列的加法與乘法為什么你的預測總在趨勢拐點翻車很多人做時間序列預測時拿到數(shù)據(jù)第一件事就是往 ARIMA 或 Prophet 里一塞跑出來一條平滑曲線看著擬合得不錯一到真實業(yè)務里就發(fā)現(xiàn)趨勢一變預測全廢。問題往往不在模型本身而在于你默認了一個假設——這個序列是「加法結構」還是「乘法結構」。這兩個詞聽起來像統(tǒng)計學課本里的老古董但在 R 語言里用 GAM廣義可加模型做時間序列分解和預測時它直接決定了你的模型是「能跟著趨勢走」還是「一到旺季就低估」。加法模型假設各成分獨立疊加觀測值 趨勢 季節(jié) 殘差。乘法模型假設成分之間是比例關系觀測值 趨勢 × 季節(jié) × 殘差。現(xiàn)實業(yè)務里銷售額隨趨勢增長時季節(jié)性波動幅度也在放大這就是典型的乘法結構。如果你硬用加法模型去擬合殘差會呈現(xiàn)明顯的異方差預測區(qū)間在高峰期窄得離譜低谷期又寬得沒用。R 語言里的 GAM 通過mgcv包提供了靈活的平滑項可以讓你在同一個框架下處理這兩種結構甚至混合結構。這篇文章面向的是已經(jīng)會用 R 做基本數(shù)據(jù)分析、但想在時間序列上把 GAM 用對的從業(yè)者。我會把加法與乘法的選擇邏輯、GAM 的平滑項設置、以及實際跑通一套分解加預測的流程講清楚順帶把那些讓我翻過車的參數(shù)坑標出來。2. 加法與乘法結構的判定先看殘差再談模型2.1 從業(yè)務場景反推結構類型拿到一條時間序列不要急著畫圖。先問一個業(yè)務問題當趨勢上升 10% 時季節(jié)波動的絕對幅度是保持不變還是也跟著漲 10%如果是前者加法結構更合理如果是后者乘法結構更合理。這個判斷不需要任何統(tǒng)計檢驗靠業(yè)務常識就能定調。舉個例子某電商平臺的日訂單量。過去兩年整體趨勢從每天 1000 單漲到 3000 單同時「雙十一」當天的峰值從 5000 單漲到 15000 單。峰值與均值的比例大致穩(wěn)定在 5 倍左右這就是乘法結構的典型信號。反過來如果是一個成熟產(chǎn)品的日活用戶趨勢基本平緩周末比工作日固定少 2000 人那加法結構就夠了。在 R 里你可以用forecast包的decompose()函數(shù)分別跑加法和乘法分解然后對比殘差圖。乘法分解的殘差如果比加法分解更接近白噪聲且殘差絕對值不隨趨勢增大而膨脹那就選乘法。library(forecast) # 假設 ts_data 是一個 ts 對象頻率為 7周季節(jié) # 加法分解 decomp_add - decompose(ts_data, type additive) # 乘法分解 decomp_mult - decompose(ts_data, type multiplicative) # 對比殘差的標準差隨趨勢的變化 # 如果乘法分解的殘差更穩(wěn)定選乘法 sd_add - sd(decomp_add$random, na.rm TRUE) sd_mult - sd(decomp_mult$random, na.rm TRUE) # 更直觀的方法看殘差與趨勢的相關系數(shù) trend_add - decomp_add$trend resid_add - decomp_add$random cor_add - cor(trend_add, abs(resid_add), use complete.obs) trend_mult - decomp_mult$trend resid_mult - decomp_mult$random cor_mult - cor(trend_mult, abs(resid_mult), use complete.obs) cat(加法殘差與趨勢相關系數(shù):, cor_add, \n) cat(乘法殘差與趨勢相關系數(shù):, cor_mult, \n)這段代碼的邏輯是如果殘差的絕對值和趨勢有正相關說明波動幅度隨趨勢增長加法模型就不合適。參數(shù)上type指定分解類型frequency在ts()里設定。注意decompose()只適合季節(jié)性固定的序列如果季節(jié)模式在變得用stl()或 GAM。2.2 用 GAM 的平滑項同時捕捉趨勢與季節(jié)GAM 的優(yōu)勢在于它不預設趨勢是線性的也不預設季節(jié)是固定的。mgcv包里的gam()函數(shù)可以用s()指定平滑項用te()指定交互項。對于時間序列我一般這樣建library(mgcv) # 構造時間索引和季節(jié)索引 n - length(ts_data) time_idx - 1:n season_idx - cycle(ts_data) # 1 到 frequency # 加法結構 GAM gam_add - gam(ts_data ~ s(time_idx, k 20) s(season_idx, k 7, bs cc), family gaussian()) # 乘法結構先取對數(shù)再跑加法 GAM gam_mult - gam(log(ts_data) ~ s(time_idx, k 20) s(season_idx, k 7, bs cc), family gaussian())這里的關鍵參數(shù)是k它控制平滑項的最大自由度。k太小會欠擬合趨勢被抹平k太大會過擬合把噪聲當信號。我一般從 10 到 30 之間試用gam.check()看殘差和k的顯著性。bs cc是循環(huán)平滑專門用于季節(jié)索引這種周期性變量不加這個12 月和 1 月之間會出現(xiàn)斷點。乘法結構取對數(shù)后跑加法 GAM等價于在原始尺度上做乘法分解。預測時記得exp()回來但要注意偏差修正——直接exp()會低估均值因為對數(shù)正態(tài)分布的均值是exp(mu sigma^2/2)。這個坑我后面會細說。2.3 模型選擇的量化依據(jù)AIC 與殘差診斷加法 GAM 和乘法 GAM 跑完后不能只看圖。AIC()可以比較但注意乘法模型是在對數(shù)尺度上算的AIC 不能直接和加法模型比。正確做法是把兩個模型的預測值都變換回原始尺度算 RMSE 或 MAE用交叉驗證比較。# 樣本外交叉驗證留出最后 30 個點 train_n - n - 30 train_data - ts_data[1:train_n] test_data - ts_data[(train_n 1):n] # 重新擬合加法模型 gam_add_train - gam(train_data ~ s(1:train_n, k 20) s(cycle(train_data), k 7, bs cc)) # 重新擬合乘法模型 gam_mult_train - gam(log(train_data) ~ s(1:train_n, k 20) s(cycle(train_data), k 7, bs cc)) # 預測 pred_add - predict(gam_add_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(test_data) )) pred_mult_log - predict(gam_mult_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(test_data) )) pred_mult - exp(pred_mult_log) # 計算 RMSE rmse_add - sqrt(mean((test_data - pred_add)^2)) rmse_mult - sqrt(mean((test_data - pred_mult)^2)) cat(加法 RMSE:, rmse_add, \n) cat(乘法 RMSE:, rmse_mult, \n)這段代碼的核心是「用樣本外誤差說話」。參數(shù)上k在訓練集上重新選不要用全量數(shù)據(jù)調好的k直接套。cycle()提取季節(jié)位置predict()的newdata必須包含和訓練時同名的變量。如果乘法模型的 RMSE 明顯更低且殘差沒有異方差那就選乘法。3. 在 R 里跑通 GAM 時間序列從數(shù)據(jù)到預測的完整鏈路3.1 數(shù)據(jù)準備與 ts 對象構造R 里做時間序列第一步是把數(shù)據(jù)轉成ts對象。很多人從 CSV 讀進來直接跑結果cycle()返回 NULL季節(jié)平滑項直接報錯。正確做法是明確頻率和起始時間。# 假設 raw_data 是數(shù)據(jù)框date 列是日期value 列是觀測值 raw_data$date - as.Date(raw_data$date) # 按日期排序 raw_data - raw_data[order(raw_data$date), ] # 構造 ts 對象頻率為 7周數(shù)據(jù) ts_data - ts(raw_data$value, frequency 7, start c(year(min(raw_data$date)), as.numeric(format(min(raw_data$date), %j)))) # 檢查 cat(頻率:, frequency(ts_data), \n) cat(周期數(shù):, length(ts_data) / frequency(ts_data), \n) head(cycle(ts_data))frequency 7表示每周 7 個觀測start參數(shù)指定起始年份和年內第幾天。如果頻率設錯比如日數(shù)據(jù)設成 365cycle()會返回 1 到 365季節(jié)平滑項k就得設得很大計算量爆炸且容易過擬合。常見做法是日數(shù)據(jù)用frequency 7捕捉周內模式或者frequency 365.25捕捉年內模式但后者需要k至少 50 以上我一般先用周頻率跑通再考慮年頻率。3.2 GAM 平滑項的參數(shù)設置與 gam.check 解讀gam()函數(shù)里最關鍵的三個參數(shù)k、bs、family。k是基函數(shù)維度決定平滑曲線的最大彎曲次數(shù)。bs是基函數(shù)類型tp是薄板回歸樣條默認cc是循環(huán)三次樣條用于周期變量cr是三次回歸樣條計算更快。family指定分布族高斯用于連續(xù)值泊松用于計數(shù)負二項用于過離散計數(shù)。# 完整模型擬合 gam_fit - gam(ts_data ~ s(time_idx, k 25, bs tp) s(season_idx, k 7, bs cc), family gaussian(), method REML) # 用 REML 估計平滑參數(shù) # 檢查 gam.check(gam_fit)gam.check()輸出四張圖殘差 vs 擬合值、QQ 圖、直方圖、響應 vs 擬合值。重點看第一張如果殘差呈現(xiàn)漏斗形隨擬合值增大而擴散說明方差非恒定加法高斯模型不合適要么換乘法取對數(shù)要么換family Gamma。QQ 圖如果尾部偏離說明殘差非正態(tài)預測區(qū)間會不準。method REML是我強烈建議的默認的 GCV 在樣本量小的時候容易過擬合REML 更穩(wěn)健。這個參數(shù)不設gam()會用 GCV殘差診斷經(jīng)常顯示k不夠其實換了 REML 就好了。3.3 預測與置信區(qū)間對數(shù)變換后的偏差修正乘法模型預測時predict()返回的是對數(shù)尺度上的值exp()回去之后得到的是中位數(shù)不是均值。如果業(yè)務要的是均值預測必須加偏差修正。# 乘法模型預測對數(shù)尺度 pred_log - predict(gam_mult, newdata data.frame( time_idx (n 1):(n 30), season_idx rep(1:7, length.out 30) ), se.fit TRUE) # 偏差修正exp(mu sigma^2/2) sigma2 - summary(gam_mult)$scale pred_mean - exp(pred_log$fit sigma2 / 2) # 置信區(qū)間近似 lower - exp(pred_log$fit - 1.96 * pred_log$se.fit) upper - exp(pred_log$fit 1.96 * pred_log$se.fit) # 注意這個區(qū)間是條件均值區(qū)間不是預測區(qū)間 # 預測區(qū)間還要加殘差方差 pred_lower - exp(pred_log$fit - 1.96 * sqrt(pred_log$se.fit^2 sigma2)) pred_upper - exp(pred_log$fit 1.96 * sqrt(pred_log$se.fit^2 sigma2))summary(gam_mult)$scale提取的是殘差方差估計。偏差修正項sigma2 / 2看起來小但當sigma2大的時候不修正會低估均值 5% 到 10%在庫存補貨場景里就是真金白銀的差距。置信區(qū)間和預測區(qū)間是兩回事前者是均值的區(qū)間后者是單次觀測的區(qū)間。業(yè)務上要「明天可能賣多少」用預測區(qū)間要「明天平均賣多少」用置信區(qū)間。4. 避坑與排查那些讓我重跑模型的參數(shù)陷阱4.1 現(xiàn)象gam.check 顯示 k 不夠調大后過擬合原因k的默認值是 10對于趨勢變化劇烈的序列10 個基函數(shù)不夠彎曲。但直接調到 50模型會把噪聲也擬合進去樣本外誤差反而上升。解決先用gam.check()看殘差是否還有模式。如果殘差已經(jīng)像白噪聲k就夠了不用管 p 值。如果殘差還有趨勢每次加 5直到殘差無模式。同時用method REML它比 GCV 更不容易過擬合。我一般從k 20開始趨勢項和季節(jié)項分開調。4.2 現(xiàn)象乘法模型預測值在低谷期為負原因對數(shù)變換后預測再exp()回去理論上不會負。但如果用了family gaussian()且沒有取對數(shù)直接跑乘法結構比如把季節(jié)項寫成乘積形式線性預測子可能為負。解決乘法結構必須在對數(shù)尺度上建模。log(ts_data)后跑加法 GAMexp()回來。如果原始數(shù)據(jù)有零或負值先加一個常數(shù)再取對數(shù)或者改用family Gamma(link log)后者直接在對數(shù)鏈接函數(shù)上建模不需要手動變換。4.3 現(xiàn)象季節(jié)平滑項 bs cc 報錯「A term has fewer unique covariate combinations than specified maximum degrees of freedom」原因k設得比季節(jié)周期的唯一值數(shù)量還大。比如周頻率frequency 7season_idx只有 7 個唯一值k設成 10 就報錯。解決k必須小于等于唯一值數(shù)量。周頻率用k 7月頻率用k 12日頻率用k 7周內模式或k 365年內模式但計算量大。如果確實需要更大的k改用bs tp并手動構造傅里葉項。4.4 現(xiàn)象預測區(qū)間在趨勢上升段明顯偏窄原因GAM 的predict(se.fit TRUE)只給了平滑項的不確定性沒有包含殘差方差。而且如果模型是加法結構但數(shù)據(jù)實際是乘法結構殘差方差會隨趨勢增大區(qū)間估計系統(tǒng)性偏窄。解決預測區(qū)間用sqrt(se.fit^2 sigma2)不要只用se.fit。同時檢查殘差是否異方差如果是換乘法結構或family Gamma。我習慣在預測后畫殘差 vs 擬合值圖確認沒有漏斗形才發(fā)結果。4.5 現(xiàn)象時間索引 time_idx 從 1 到 n預測時新數(shù)據(jù)的時間索引對不上原因訓練時time_idx是 1 到 n預測時新數(shù)據(jù)的time_idx必須從 n1 開始。如果直接用1:30模型會以為你在預測歷史區(qū)間。解決預測時構造newdata的time_idx用(n1):(nh)season_idx用cycle()或手動指定。更穩(wěn)妥的做法是把時間索引存成數(shù)據(jù)框的一列訓練和預測用同一套構造邏輯。5. 進階技巧用 te() 捕捉趨勢與季節(jié)的交互效應5.1 什么時候需要交互項加法 GAM 假設趨勢和季節(jié)是獨立作用的趨勢上升不影響季節(jié)波動的形狀。但現(xiàn)實中很多序列的季節(jié)模式會隨趨勢變化。比如一個產(chǎn)品剛上市時只有周末賣得好隨著市場成熟工作日銷量也上來了季節(jié)波動的幅度和形狀都變了。這時候s(time_idx) s(season_idx)就不夠了需要用te(time_idx, season_idx)張量積平滑來捕捉交互。# 張量積交互模型 gam_te - gam(ts_data ~ te(time_idx, season_idx, k c(15, 7), bs c(tp, cc)), family gaussian(), method REML) # 對比無交互模型 gam_no_te - gam(ts_data ~ s(time_idx, k 15) s(season_idx, k 7, bs cc), family gaussian(), method REML) # 用 AIC 比較同一尺度下 cat(交互模型 AIC:, AIC(gam_te), \n) cat(無交互 AIC:, AIC(gam_no_te), \n)te()的參數(shù)k是一個向量分別對應兩個變量的基函數(shù)維度。bs也是向量對應各自的基函數(shù)類型。交互模型的自由度是k[1] * k[2]的量級計算量比加法模型大很多樣本量少于 200 時慎用容易過擬合。我一般先跑加法模型如果殘差在特定時間段比如旺季還有系統(tǒng)偏差再上交互項。5.2 用 vis.gam 可視化交互效應mgcv自帶的vis.gam()可以畫三維透視圖或等高線圖直觀看到季節(jié)模式如何隨趨勢變化。# 三維透視圖 vis.gam(gam_te, view c(time_idx, season_idx), theta 30, phi 20, color heat, plot.type persp, ticktype detailed) # 等高線圖 vis.gam(gam_te, view c(time_idx, season_idx), plot.type contour, color heat)view指定兩個變量theta和phi控制視角plot.type選persp或contour。等高線圖里如果等高線是平行的說明沒有交互如果彎曲或交叉說明季節(jié)模式隨趨勢變了。這個圖我每次做完交互模型都會看比 AIC 更直觀。5.3 一個我常犯的錯誤交互項加了但沒做預測對比剛用te()的時候我看到 AIC 降了就以為交互模型更好直接拿去預測結果樣本外 RMSE 反而更高。后來養(yǎng)成習慣不管 AIC 多低一定做樣本外交叉驗證。交互模型的方差更大樣本量不夠時AIC 的懲罰項不足以抵消過擬合。# 樣本外對比 train_n - n - 30 gam_te_train - gam(ts_data[1:train_n] ~ te(1:train_n, cycle(ts_data[1:train_n]), k c(15, 7), bs c(tp, cc)), family gaussian(), method REML) pred_te - predict(gam_te_train, newdata data.frame( time_idx (train_n 1):n, season_idx cycle(ts_data[(train_n 1):n]) )) rmse_te - sqrt(mean((ts_data[(train_n 1):n] - pred_te)^2)) cat(交互模型樣本外 RMSE:, rmse_te, \n)如果rmse_te比加法模型的樣本外 RMSE 還大說明交互項在擬合噪聲果斷回退。這個習慣讓我少了很多「AIC 好看但業(yè)務翻車」的情況。5.4 最后的習慣先畫圖再跑模型我現(xiàn)在拿到任何時間序列第一件事是plot(ts_data)第二件事是plot(decompose(ts_data))或plot(stl(ts_data, s.window periodic))。圖上看一眼趨勢和季節(jié)的形態(tài)比任何統(tǒng)計檢驗都快。加法還是乘法很多時候圖上一眼就能定如果季節(jié)波動的幅度隨趨勢明顯放大乘法如果幅度穩(wěn)定加法。GAM 只是把這個判斷量化并給出預測區(qū)間。希望幫到你。本文還有配套的精品資源點擊獲取