Bagaimana menemukan model semi-sinusoidal yang cocok untuk R?

Saya ingin berasumsi bahwa suhu permukaan laut dari Laut Baltik adalah tahun yang sama tahun demi tahun, dan kemudian menggambarkannya dengan model fungsi / linier. Gagasan yang saya miliki adalah hanya memasukkan tahun sebagai angka desimal (atau num_months / 12) dan keluar berapa suhu yang seharusnya sekitar waktu itu. Melemparkannya ke fungsi lm () di R, itu tidak mengenali data sinusoidal sehingga hanya menghasilkan garis lurus. Jadi saya meletakkan fungsi sin () dalam braket I () dan mencoba beberapa nilai agar sesuai dengan fungsi tersebut, dan itu mendekati apa yang saya inginkan. Tetapi laut semakin cepat memanas di musim panas dan kemudian mendingin lebih lambat di musim gugur ... Jadi modelnya salah pada tahun pertama, kemudian menjadi lebih benar setelah beberapa tahun, dan kemudian di masa depan saya kira itu menjadi lebih dan lebih salah lagi.

Bagaimana saya bisa mendapatkan R untuk memperkirakan model untuk saya, jadi saya tidak perlu menebak angka sendiri? Kuncinya di sini adalah bahwa saya ingin menghasilkan nilai yang sama tahun demi tahun, tidak hanya benar untuk satu tahun. Jika saya tahu lebih banyak tentang matematika, mungkin saya bisa menebaknya sebagai sesuatu seperti Poisson atau Gaussian alih-alih dosa (), tapi saya juga tidak tahu bagaimana melakukannya. Setiap bantuan untuk mendekati jawaban yang baik akan sangat dihargai.

Berikut adalah data yang saya gunakan, dan kode untuk menunjukkan hasil sejauh ini:

# SST from Bradtke et al 2010
ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)
SSTlm <- lm(SST$Degrees ~ I(sin(pi*2.07*SST$ToY)))
summary(SSTlm)
plot(SST,xlim=c(0,4),ylim=c(0,17))
par(new=T)
plot(data.frame(ToY=SST$ToY,Degrees=8.4418-6.9431*sin(2.07*pi*SST$ToY)),type="l",xlim=c(0,4),ylim=c(0,17))

r regression time-series lm GaRyu
sumber

Jawaban:

Ini dapat dilakukan dengan regresi linier -

Anda hanya perlu baik dan istilah di setiap frekuensi. $\sin$ $\cos$

Alasan mengapa Anda dapat menggunakan istilah dan dalam regresi linier untuk menangani musiman dengan amplitudo dan fase apa pun adalah karena identitas trigonometri berikut : $\sin$ $\cos$

Gelombang sinus 'umum' dengan dengan amplitudo dan fase , , dapat ditulis sebagai kombinasi linear mana dan sedemikian rupa sehingga $A$ $\varphi$ $A \sin (x + \varphi)$ $a\sin x+b\cos x$ $a$ $b$ dan $A=\sqrt{a^2+b^2}$ . Mari kita lihat bahwa keduanya setara: $\sin\varphi = \frac{b}{\sqrt{a^2+b^2}}$

\begin{array}{rcl} Sebuah dosa (x) + b \cos (x) & = & \sqrt{{Sebuah}^{2} + b^{2}} (\frac{Sebuah}{\sqrt{{Sebuah}^{2} + b^{2}}} dosa (x) + \frac{b}{\sqrt{{Sebuah}^{2} + b^{2}}} \cos (x)) \\ = & SEBUAH [dosa (x) \cos (φ) + \cos (x) dosa (φ)] \\ = & SEBUAH dosa (x + φ) . \end{array}

$\begin{eqnarray} a \sin(x) + b \cos(x) &=& \sqrt{a^2+b^2} \left(\frac{a}{\sqrt{a^2+b^2}} \sin(x) + \frac{b}{\sqrt{a^2+b^2}} \cos(x)\right)\\ &=& A\left[\sin(x)\cos(\varphi) + \cos(x)\sin(\varphi)\right]\\ &=& A\sin(x+\varphi)\,\text{.} \end{eqnarray}$

Inilah model 'dasar':

 SSTlm <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY),data=SST)
 summary(SSTlm)

[menggunting]

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)              8.292      0.135   61.41   <2e-16 *** 
sin(2 * pi * ToY)       -5.916      0.191  -30.98   <2e-16 ***  
cos(2 * pi * ToY)       -4.046      0.191  -21.19   <2e-16 *** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9355 on 45 degrees of freedom
Multiple R-squared: 0.969,      Adjusted R-squared: 0.9677 
F-statistic: 704.3 on 2 and 45 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylim=c(1.5,16.5),data=SST)
 lines(SST$ToY,SSTlm$fitted,col=2)

dosa cocok

Sunting: Catatan penting - Istilah bekerja karena periode fungsi telah diatur sehingga satu periode = 1 unit . Jika periode berbeda dari 1, misalkan periodenya , maka Anda perlu $2\pi\,t$ $t$ $\omega$ sebagai gantinya. $(2\pi/\omega)\, t$

Inilah model dengan harmonik kedua:

 SSTlm2 <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY)
                        +sin(4*pi*ToY)+cos(4*pi*ToY),data=SST)
 summary(SSTlm2)

[menggunting]

Coefficients:
                  Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        8.29167    0.02637  314.450  < 2e-16 ***  
sin(2 * pi * ToY) -5.91562    0.03729 -158.634  < 2e-16 ***  
cos(2 * pi * ToY) -4.04632    0.03729 -108.506  < 2e-16 ***  
sin(4 * pi * ToY)  1.21244    0.03729   32.513  < 2e-16 ***  
cos(4 * pi * ToY)  0.33333    0.03729    8.939 2.32e-11 ***  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1827 on 43 degrees of freedom
Multiple R-squared: 0.9989,     Adjusted R-squared: 0.9988 
F-statistic:  9519 on 4 and 43 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylab="Degrees",xlab="ToY",ylim=c(1.5,16.5),data=SST)
 lines(SSTlm2$fitted~ToY,col=2,data=SST)

dosa cocok 2

... dan seterusnya, dengan 6*pi*ToYdll. Jika ada sedikit kebisingan dalam data saya mungkin akan berhenti dengan model kedua ini.

Dengan persyaratan yang cukup, Anda dapat dengan tepat menyesuaikan urutan periodik yang asimetris dan bahkan bergerigi, tetapi kecocokan yang dihasilkan mungkin 'bergoyang'. Berikut adalah fungsi asimetris (ini gigi gergaji - ) ditambahkan ke versi skala fungsi periodik Anda), dengan harmonik ketiga (merah) dan keempat (hijau). Fit hijau rata-rata sedikit lebih dekat tetapi "goyah" (bahkan ketika fit melewati setiap titik, cocok mungkin sangat goyah di antara poin).

sin fit 3 & 4

$\cos$ $\sin$

Jika Anda ingin pas yang lebih halus daripada yang dihasilkan oleh pendekatan ini pada seri yang tidak mulus, Anda mungkin ingin melihat cocok spline periodik .

Namun pendekatan lain adalah dengan menggunakan boneka musiman, tetapi pendekatan sin / cos sering lebih baik jika itu adalah fungsi periodik yang lancar.

Pendekatan semacam ini untuk musiman juga dapat beradaptasi dengan situasi di mana musiman berubah, seperti menggunakan trigonometrik atau musiman boneka dengan model ruang negara.

Sementara pendekatan model linier yang dibahas di sini mudah digunakan, satu keuntungan dari pendekatan regresi nonlinear @ COOLSerdash adalah bahwa ia dapat menangani berbagai situasi yang jauh lebih luas - Anda tidak perlu banyak berubah sebelum Anda berada dalam situasi di mana linier regresi tidak lagi cocok tetapi kuadrat terkecil nonlinear masih dapat digunakan (memiliki periode yang tidak diketahui akan menjadi salah satu kasus seperti itu).

Glen_b -Reinstate Monica
sumber

Luar biasa! Terima kasih, saya benar-benar harus mencoba mempelajari lebih lanjut tentang metode untuk menangani frekuensi. Saya tidak begitu mengerti mengapa bagian cos diperlukan, tetapi mengetahui prinsipnya membuatnya mudah untuk diterapkan.

GaRyu

@COOLSerdash - sebenarnya, saya harap Anda belum menghapus jawaban Anda (memang saya telah memutarnya); itu memiliki keuntungan bekerja dalam berbagai keadaan yang jauh lebih luas; tweak beberapa hal tentang masalah dan Anda bisa kehilangan linearitas - dan kemudian pendekatan saya tidak berguna, tetapi Anda masih berfungsi. Saya pikir ada banyak yang bisa dikatakan untuk bisa melakukannya dengan cara itu.

Glen_b -Reinstate Monica

@ Glen_b Ah maaf, saya pikir posting Anda membuat saya menjadi berlebihan karena saya tidak menggunakan cara standar untuk menangani masalah tersebut. Saya menghapusnya.

COOLSerdash

\cos

$\cos$

Itu bukan saya .... Anda mengatakan fase offset seolah-olah itu bernama apa yang sedang terjadi, dan itu secara matematis. Tetapi bagi Anda poin kuncinya adalah bahwa kemungkinan 31 Desember / 1 Januari adalah asal sewenang-wenang untuk sepanjang tahun karena keterlambatan dalam respons suhu terhadap variasi dalam penerimaan radiasi. Jadi fase offset adalah nama di sini untuk sesuatu klimatologis juga, waktu suhu minimum dan maksimum relatif terhadap sistem perekaman Anda. (Ini detail kecil tapi saya lebih suka menghitung waktu tahun selama 12 bulan sebagai 1/24, 3/24, ..., 23/24.)

Nick Cox

Suhu yang Anda berikan dalam pertanyaan Anda berulang setiap tahun. Saya menduga ini bukan suhu yang benar-benar diukur selama empat tahun. Dalam contoh Anda, Anda tidak perlu model, karena suhunya hanya mengulangi persis. Tetapi sebaliknya Anda bisa menggunakan nlsfungsi ini agar sesuai dengan kurva sinus:

ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)

par(cex=1.5, bg="white")
plot(Degrees~ToY,xlim=c(0,4),ylim=c(0,17), pch=16, las=1)

nls.mod <-nls(Degrees ~ a + b*sin(2*pi*c*ToY), start=list(a = 1, b = 1, c=1))

co <- coef(nls.mod) 
f <- function(x, a, b, c) {a + b*sin(2*pi*c*x) }

curve(f(x, a=co["a"], b=co["b"], c=co["c"]), add=TRUE ,lwd=2, col="steelblue")

NLS cocok

Tapi pas tidak terlalu bagus, terutama di awal. Tampaknya data Anda tidak dapat dimodelkan secara memadai oleh kurva sinus sederhana. Mungkin fungsi trigonometri yang lebih kompleks akan berhasil?

nls.mod2 <-nls(Degrees ~ a + b*sin(2*pi*c*ToY)+d*cos(2*pi*e*ToY), start=list(a = 1, b = 1, c=1, d=1, e=1))

co2 <- coef(nls.mod2) 
f <- function(x, a, b, c, d, e) {a + b*sin(2*pi*c*x)+d*cos(2*pi*e*x) }

curve(f(x, a=co2["a"], b=co2["b"], c=co2["c"], d=co2["d"], e=co2["e"]), add=TRUE ,lwd=2, col="red")

NLS cocok 2

Kurva merah lebih cocok dengan data. Dengan nlsfungsinya, Anda bisa memasukkan model yang menurut Anda sesuai.

Atau mungkin Anda bisa menggunakan forecastpaket itu. Dalam contoh di bawah ini, saya berasumsi bahwa rangkaian waktu dimulai pada Januari 2010:

library(forecast)

Degrees.ts <- ts(Degrees, start=c(2010,1), frequency=12)

Degree.trend <- auto.arima(Degrees.ts)

degrees.forecast <- forecast(Degree.trend, h=12, level=c(80,95), fan=F)

plot(degrees.forecast, las=1, main="", xlab="Time", ylab="Degrees")

ARIMA

Karena data bersifat deterministik, tidak ada pita kepercayaan yang ditampilkan.

COOLSerdash
sumber

Tidak ada alasan untuk kuadrat-nonlinier di sini, bukan karena itu tidak akan berfungsi dengan baik. Hitung sin (2 * pi * TOY), cos (2 * pi * ToY) di muka dan beri mereka makan lm()seperti prediktor lain. Dengan kata lain, lm()tidak perlu melihat trigonometri sama sekali. Namun, Anda mungkin perlu model lain untuk menangkap asimetri dengan baik. Saya bukan pengguna R biasa tapi saya sering menggunakan pendekatan ini di tempat lain (lihat stata-journal.com/sjpdf.html?articlenum=st0116 ).

Nick Cox

@NickCox Terima kasih Nick, itu saran yang sangat membantu. Saya akan memperbarui jawaban saya sedikit.

COOLSerdash

Glen lebih cepat :)

COOLSerdash

@COOLserdash Saya bahkan tidak melihat komentar Nick Cox di sana; itu datang ketika saya menghasilkan jawaban saya. (Pendekatan ini cukup jelas jika Anda pernah melihat seri Fourier.)

Glen_b -Reinstate Monica

Seperti yang disiratkan @Glen_b, ini adalah pendekatan standar, tidak diketahui secara universal.

Nick Cox