Sesuaikan istilah sinusoidal dengan data

Meskipun saya membaca posting ini , saya masih tidak tahu bagaimana menerapkan ini pada data saya sendiri dan berharap seseorang dapat membantu saya.

Saya memiliki data berikut:

y <- c(11.622967, 12.006081, 11.760928, 12.246830, 12.052126, 12.346154, 12.039262, 12.362163, 12.009269, 11.260743, 10.950483, 10.522091,  9.346292,  7.014578,  6.981853,  7.197708,  7.035624,  6.785289, 7.134426,  8.338514,  8.723832, 10.276473, 10.602792, 11.031908, 11.364901, 11.687638, 11.947783, 12.228909, 11.918379, 12.343574, 12.046851, 12.316508, 12.147746, 12.136446, 11.744371,  8.317413, 8.790837, 10.139807,  7.019035,  7.541484,  7.199672,  9.090377,  7.532161,  8.156842,  9.329572, 9.991522, 10.036448, 10.797905)
t <- 18:65

Dan sekarang saya hanya ingin menyesuaikan gelombang sinus

dengan empat tidak diketahui , , dan untuk itu.ω ϕ $A$$\omega$$\phi$$C$

Sisa dari kode saya terlihat adalah sebagai berikut

res <- nls(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=list(A=1,omega=1,phi=1,C=1))
co <- coef(res)

fit <- function(x, a, b, c, d) {a*sin(b*x+c)+d}

# Plot result
plot(x=t, y=y)
curve(fit(x, a=co["A"], b=co["omega"], c=co["phi"], d=co["C"]), add=TRUE ,lwd=2, col="steelblue")

Tetapi hasilnya benar-benar buruk.

Saya akan sangat menghargai bantuan apa pun.

Tepuk tangan.

Jika Anda hanya ingin perkiraan yang baik dari dan tidak terlalu peduli dengan kesalahan standarnya:$\omega$

ssp <- spectrum(y)  
per <- 1/ssp$freq[ssp$spec==max(ssp$spec)]
reslm <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t))
summary(reslm)

rg <- diff(range(y))
plot(y~t,ylim=c(min(y)-0.1*rg,max(y)+0.1*rg))
lines(fitted(reslm)~t,col=4,lty=2)   # dashed blue line is sin fit

# including 2nd harmonic really improves the fit
reslm2 <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t))
summary(reslm2)
lines(fitted(reslm2)~t,col=3)    # solid green line is periodic with second harmonic

(Kecocokan yang lebih baik mungkin akan menjelaskan pencilan dalam seri itu dalam beberapa cara, mengurangi pengaruh mereka.)

---

Jika Anda menginginkan gagasan tentang ketidakpastian di , Anda dapat menggunakan kemungkinan profil ( pdf1 , pdf2 - referensi untuk mendapatkan perkiraan CI atau UK dari kemungkinan profil atau variannya tidak sulit ditemukan)$\omega$

(Atau, Anda dapat memberi makan perkiraan ini ke nls ... dan memulainya yang sudah konvergen.)

$2 \pi / 20$

Ketika saya memasukkannya ke nlsdalam startdaftar, saya mendapatkan kurva yang jauh lebih masuk akal, meskipun masih memiliki beberapa bias sistematis.

Bergantung pada apa tujuan Anda dengan kumpulan data ini, Anda dapat mencoba meningkatkan kecocokan dengan menambahkan istilah tambahan atau menggunakan pendekatan nonparametrik seperti proses Gaussian dengan kernel berkala.

Memilih nilai awal secara otomatis

Jika Anda ingin memilih frekuensi yang dominan, Anda dapat menggunakan fast Fourier transform (FFT). Ini adalah jalan keluar dari bidang keahlian saya, jadi saya akan membiarkan orang lain mengisi rincian jika mereka mau (terutama tentang langkah 2 dan 3), tetapi Rkode di bawah ini akan berfungsi.

# Step 1: do the FFT
raw.fft = fft(y)

# Step 2: drop anything past the N/2 - 1th element.
# This has something to do with the Nyquist-shannon limit, I believe
# (https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem)
truncated.fft = raw.fft[seq(1, length(y)/2 - 1)]

# Step 3: drop the first element. It doesn't contain frequency information.
truncated.fft[1] = 0

# Step 4: the importance of each frequency corresponds to the absolute value of the FFT.
# The 2, pi, and length(y) ensure that omega is on the correct scale relative to t.
# Here, I set omega based on the largest value using which.max().
omega = which.max(abs(truncated.fft)) * 2 * pi / length(y)

Anda juga dapat merencanakan abs(truncated.fft)untuk melihat apakah ada frekuensi penting lainnya, tetapi Anda harus sedikit bermain-main dengan penskalaan sumbu-x.

Juga, saya percaya @Glen_b benar bahwa masalahnya cembung setelah Anda tahu omega (atau mungkin Anda perlu tahu phi juga? Saya tidak yakin). Bagaimanapun, mengetahui nilai awal untuk parameter lain seharusnya tidak sepenting untuk omega jika mereka berada di stadion baseball yang tepat. Anda mungkin bisa mendapatkan perkiraan yang layak dari parameter lain dari FFT, tapi saya tidak yakin bagaimana cara kerjanya.

Sebagai alternatif dari apa yang telah dikatakan, mungkin perlu dicatat bahwa model AR (2) dari kelas model ARIMA dapat digunakan untuk menghasilkan perkiraan dengan pola gelombang sinus.

Panratz (1991) memberi tahu kita hal berikut tentang siklus stokastik:

Pola siklus stokastik dapat dianggap sebagai pola gelombang sinus terdistorsi dalam pola perkiraan: Ini adalah gelombang sinus dengan periode stokastik (probabilistik), amplitudo, dan sudut fase.

Untuk melihat apakah model seperti itu dapat dipasang ke data saya menggunakan auto.arima()fungsi dari paket perkiraan untuk mengetahui apakah itu akan menyarankan model AR (2). Ternyata auto.arima()fungsi menyarankan model ARMA (2,2); bukan model AR (2) murni, tapi ini OK. Tidak apa-apa karena model ARMA (2,2) berisi komponen AR (2), jadi aturan yang sama (tentang siklus stokastik) berlaku. Artinya, kita masih dapat memeriksa kondisi tersebut untuk melihat apakah prakiraan gelombang sinus akan diproduksi.

Hasilnya auto.arima(y)ditunjukkan di bawah ini.

Series: y 
ARIMA(2,0,2) with non-zero mean 

Coefficients:
         ar1      ar2      ma1     ma2  intercept
      1.7347  -0.8324  -1.2474  0.6918    10.2727
s.e.  0.1078   0.0981   0.1167  0.1911     0.5324

sigma^2 estimated as 0.6756:  log likelihood=-60.14
AIC=132.27   AICc=134.32   BIC=143.5

Plot di bawah ini menunjukkan seri asli, y, kesesuaian model ARMA (2,2), dan 14 perkiraan out-of-sample. Seperti yang dapat dilihat, perkiraan out-of-sample mengikuti pola gelombang sinus.

Ingat dua hal. 1) Ini hanya analisis yang sangat cepat (menggunakan alat otomatis) dan perawatan yang tepat akan melibatkan mengikuti metodologi Box-Jenkins. 2) Perkiraan ARIMA bagus dalam peramalan jangka pendek, sehingga Anda dapat menemukan bahwa ramalan jangka panjang dari model dalam jawaban oleh @ David J. Harris dan @Glen_b agar lebih andal.

Terakhir, semoga ini adalah tambahan yang bagus untuk beberapa jawaban yang sudah sangat informatif.

Referensi : Peramalan dengan model regresi dinamis: Alan Pankratz, 1991, (John Wiley and Sons, New York), ISBN 0-471-61528-5

Metode saat ini untuk mencocokkan kurva dosa ke set data yang diberikan membutuhkan tebakan pertama parameter, diikuti oleh proses interatif. Ini adalah masalah regresi non-linear. Metode yang berbeda terdiri dari transformasi regresi non-linear menjadi regresi linier berkat persamaan integral yang nyaman. Kemudian, tidak perlu untuk menebak awal dan tidak perlu untuk proses berulang: pemasangan langsung diperoleh. Dalam hal fungsi y = a + r * sin (w * x + phi) atau y = a + b * sin (w * x) + c * cos (w * x), lihat halaman 35-36 kertas "Régress sinusoidale" diterbitkan di Scribd: http://www.scribd.com/JJacquelin/documents Dalam hal fungsi y = a + p * x + r * sin (w * x + phi): halaman 49-51 dari bab "Regresi linier dan sinusoidal campuran". Dalam hal fungsi yang lebih rumit, proses umum dijelaskan dalam bab "Regenerasi sinusoidal umum" halaman 54-61, diikuti oleh contoh numerik y = r * sin (w * x + phi) + (b / x) + c * ln (x), halaman 62-63

Jika Anda mengetahui titik terendah dan tertinggi dari data yang tampak kosinus, Anda dapat menggunakan fungsi sederhana ini untuk menghitung semua koefisien kosinus:

getMyCosine <- function(lowest_point=c(pi,-1), highest_point=c(0,1)){
  cosine <- list(
    T = pi / abs(highest_point[1] - lowest_point[1]),
    b = - highest_point[1],
    k = (highest_point[2] + lowest_point[2]) / 2,
    A = (highest_point[2] - lowest_point[2]) / 2
  )
  return(cosine)
}

Di bawahnya digunakan untuk mensimulasikan variasi suhu sepanjang hari dengan fungsi kosinus, dengan memasukkan nilai jam dan suhu untuk jam terendah dan terhangat:

c <- getMyCosine(c(4,10),c(17,25)) 
# lowest temprature at 4:00 (10 degrees), highest at 17:00 (25 degrees)

x = seq(0,23,by=1);  y = c$A*cos(c$T*(x +c$b))+c$k ; 
library(ggplot2);   qplot(x,y,geom="step")

Outputnya di bawah ini:

Pilihan lain adalah menggunakan fungsi generik optim atau nls. Saya sudah mencoba keduanya tidak ada yang benar-benar kuat

Fungsi-fungsi berikut ini mengambil data dalam y dan menghitung parameter.

calc.period <- function(y,t)
{     
   fs <- 1/(t[2]-t[1])
   ssp <- spectrum(y,plot=FALSE )  
   fN <- ssp$freq[which.max(ssp$spec)]
   per <- 1/(fN*fs)
   return(per)
 }

fit.sine<- function(y, t)
{ 
  data <- data.frame(x = as.vector(t), y=as.vector(y))
  min.RSS <- function (data, par){
    with(data, sum((par[1]*sin(2*pi*par[2]*x + par[3])+par[4]-y )^2))
  }  
  amp = sd(data$y)*2.**0.5
  offset = mean(data$y)
  fest <- 1/calc.period(y,t)
  guess = c( amp, fest,  0,   offset)
  #res <- optim(par=guess, fn = min.RSS, data=data ) 
  r<-nls(y~offset+A*sin(2*pi*f*t+phi), 
     start=list(A=amp, f=fest, phi=0, offset=offset))
  res <- list(par=as.vector(r$m$getPars()))
  return(res)
}

 genSine <- function(t, params)
     return( params[1]*sin(2*pi*params[2]*t+ params[3])+params[4])

penggunaannya adalah sebagai berikut:

t <- seq(0, 10, by = 0.01)
A <- 2 
f <- 1.5
phase <- 0.2432
offset <- -2

y <- A*sin(2*pi*f*t +phase)+offset + rnorm(length(t), mean=0, sd=0.2)

reslm1 <- fit.sine(y = y, t= t)

Kode berikut membandingkan data

ysin <- genSine(as.vector(t), params=reslm1$par)
ysin.cor <- genSine(as.vector(t), params=c(A, f, phase, offset))

plot(t, y)
lines(t, ysin, col=2)
lines(t, ysin.cor, col=3)