Passen Sie einen sinusförmigen Term an Daten an

26

Obwohl ich diesen Beitrag gelesen habe, weiß ich immer noch nicht, wie ich das auf meine eigenen Daten anwenden soll, und hoffe, dass mir jemand helfen kann.

Ich habe folgende Daten:

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

Und jetzt möchte ich einfach eine Sinuswelle anpassen

y(t)=Asin(ωt+ϕ)+C.

mit den vier Unbekannten A , ω , ϕ und C dazu.

Der Rest meines Codes sieht folgendermaßen aus

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")

Aber das Ergebnis ist wirklich schlecht.

Sinus fit

Ich würde mich über jede Hilfe sehr freuen.

Prost.

Pascal
quelle
Sie versuchen, eine Sinuswelle an die Daten anzupassen, oder Sie versuchen, eine Art harmonisches Modell mit einer Sinus- und einer Kosinuskomponente anzupassen? Es gibt eine harmonische Funktion im TSA-Paket in R, die Sie möglicherweise überprüfen möchten. Passen Sie Ihr Modell damit an und sehen Sie, welche Ergebnisse Sie erzielen.
Eric Peterson
5
Haben Sie verschiedene Startwerte ausprobiert? Ihre Verlustfunktion ist nicht konvex, sodass unterschiedliche Ausgangswerte zu unterschiedlichen Lösungen führen können.
Stefan Wager
1
Erzählen Sie uns mehr über die Daten. In der Regel gibt es eine bekannte Periodizität, so dass nicht aus den Daten abgeschätzt werden muss. Ist das eine Zeitreihe oder etwas anderes? Es ist viel einfacher, wenn Sie durch ein lineares Modell getrennte Sinus- und Cosinus-Terme anpassen können.
Nick Cox
2
Wenn Sie einen unbekannten Zeitraum haben, wird Ihr Modell nichtlinear (auf ein solches Ereignis wird in der ausgewählten Antwort im verknüpften Beitrag hingewiesen). Vorausgesetzt, die anderen Parameter sind bedingt linear; Für einige nichtlineare LS-Routinen sind diese Informationen wichtig und können das Verhalten verbessern. Eine Möglichkeit könnte darin bestehen, spektrale Methoden zu verwenden, um die Periode und die Bedingung dafür zu ermitteln. Eine andere Möglichkeit besteht darin, die Periode und die anderen Parameter über eine nichtlineare bzw. lineare Optimierung iterativ zu aktualisieren.
Glen_b -Reinstate Monica
(I bearbeiten nur die Antwort gibt den besonderen Fall des unbekannten Zeitraums machen ein explizites Beispiel dafür , was machen kann es nicht linear.)
Glen_b -Reinstate Monica

Antworten:

18

Wenn Sie nur eine gute Schätzung von wünschen und sich nicht für den Standardfehler interessieren:ω

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

Sinusplot

(Eine noch bessere Anpassung würde vielleicht die Ausreißer in dieser Reihe in irgendeiner Weise erklären und ihren Einfluss verringern.)

---

Wenn Sie eine Vorstellung von der Unsicherheit in möchten , können Sie die Profilwahrscheinlichkeit verwenden ( pdf1 , pdf2 - Referenzen zum Abrufen von ungefähren CIs oder SEs aus der Profilwahrscheinlichkeit oder ihren Varianten sind nicht schwer zu finden).ω

(Alternativ können Sie diese Schätzungen in nls eingeben ... und es bereits konvergiert starten.)

Glen_b - Setzen Sie Monica wieder ein
quelle
(+1) nette Antwort. Ich habe versucht, das lineare Modell mit anzupassen, lm(y~sin(2*pi*t)+cos(2*pi*t)aber das hat nicht funktioniert ( cosTerm war immer 1). Nur aus Neugier: Was machen die ersten beiden Linien (ich weiß, dass sie spectrumdie spektrale Dichte schätzen)?
COOLSerdash
1
t , die der Zeitraum sind (wie in der verknüpften Frage angegeben), um 2*pi*tzu arbeiten. Ich sollte zurückgehen und das in der anderen Antwort betonen. (ctd)
Glen_b -Reinstate Monica
1
@COOLSerdash (ctd) - Die 2. Zeile ermittelt die Frequenz, die dem größten Peak im Spektrum zugeordnet ist, und kehrt sie um, um die Periode zu identifizieren. Zumindest in diesem Fall (ich vermute jedoch, dass es sich um einen weiter gefassten Standard handelt) identifiziert der Standard den Zeitraum, der die Wahrscheinlichkeit so stark maximiert, dass ich die Schritte gelöscht habe, die ich zur Maximierung der Profilwahrscheinlichkeit in der Region um diesen Zeitraum durchgeführt habe. Die Funktion specin TSA ist vielleicht besser (es scheint mehr Optionen zu geben, von denen eine manchmal wichtig sein kann), aber in diesem Fall war der Hauptpeak genau an der gleichen Stelle wie bei, spectrumalso habe ich mich nicht darum gekümmert.
Glen_b
@ Glen_b diese Methode wirkt Wunder für meinen Anwendungsfall. Ich muss auch eine cos (x) -Kurve anpassen, aber es funktioniert nicht so gut ... Ich habe die auf geändert reslm, reslm <- lm(y ~ cos(2*pi/per*t)+tan(2*pi/per*t))aber das sieht nicht richtig aus. irgendwelche Hinweise?
Amit Kohli
Warum hast du dort eine Beige-Benennung?
Glen_b
15

2π/20

Als ich das auf nlsdie startListe setzte, bekam ich eine Kurve, die viel vernünftiger war, obwohl sie immer noch einige systematische Vorurteile aufweist.

Je nachdem, welches Ziel Sie mit diesem Datensatz verfolgen, können Sie versuchen, die Anpassung zu verbessern, indem Sie zusätzliche Terme hinzufügen oder einen nichtparametrischen Ansatz wie einen Gaußschen Prozess mit einem periodischen Kernel verwenden.

Sinus fit

Startwert automatisch wählen

Wenn Sie die dominante Frequenz auswählen möchten, können Sie eine schnelle Fourier-Transformation (FFT) verwenden. Dies liegt außerhalb meines Fachgebiets, sodass andere Leute die Details eintragen können, wenn sie dies wünschen (insbesondere zu den Schritten 2 und 3), aber der folgende RCode sollte funktionieren.

# 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)

Sie können auch zeichnen, um abs(truncated.fft)zu sehen, ob es andere wichtige Frequenzen gibt, aber Sie müssen ein wenig mit der Skalierung der x-Achse herumspielen.

Außerdem glaube ich, dass @Glen_b richtig ist, dass das Problem konvex ist, sobald Sie Omega kennen (oder müssen Sie vielleicht auch Phi kennen? Ich bin nicht sicher). In jedem Fall sollte die Kenntnis der Startwerte für die anderen Parameter nicht annähernd so wichtig sein wie für Omega, wenn sie sich im richtigen Ballpark befinden. Sie könnten wahrscheinlich anständige Schätzungen der anderen Parameter von der FFT erhalten, aber ich bin nicht sicher, wie das funktionieren würde.

David J. Harris
quelle
1
Danke für diesen Hinweis. Um es ein wenig zu verdeutlichen: Die Daten sind Teil eines Microarrays, bei dem die Periodizität von Genen über die Zeit gemessen wurde, dh die gezeigten Daten sind die Expressionsdaten eines Gens. Das Problem ist jetzt, dass ich diese Methode auf ungefähr 40.000 Gene anwenden möchte, die alle unterschiedliche Periodizitäten und Amplituden aufweisen. Es ist also sehr wichtig, dass unabhängig von den Ausgangsbedingungen eine gute Passform gefunden wird.
Pascal
1
@ Pascal In meinen Updates oben finden Sie eine Empfehlung für die automatische Auswahl des Startwerts für Omega.
David J. Harris
2
ϕab
Ich frage mich, wo die x-Werte hier ins Spiel kommen. Sicher macht es einen Unterschied für Omega, ob die angegebenen y-Werte durch 1 oder durch 5 x-Schritte getrennt sind, nicht wahr?
Knub
1
Programmiertipp ohne Bezug zur Frage: Vorsicht bei der Benennung von R-Objekten als foo.bar. Dies liegt daran, wie R Methoden für Klassen angibt .
Firebug
10

Als Alternative zu dem, was bereits gesagt wurde, kann angemerkt werden, dass ein AR (2) -Modell aus der Klasse der ARIMA-Modelle verwendet werden kann, um Vorhersagen mit einem Sinuswellenmuster zu generieren.

yt=C+ϕ1yt1+ϕ2yt2+at
Cϕ1ϕ2at

ϕ12+4ϕ2<0.

Panratz (1991) erzählt uns folgendes über stochastische Zyklen:

Ein stochastisches Zyklusmuster kann als verzerrtes Sinuswellenmuster im Prognosemuster betrachtet werden: Es ist eine Sinuswelle mit einer stochastischen (probabilistischen) Periode, Amplitude und einem Phasenwinkel.

Um zu sehen, ob ein solches Modell an die Daten angepasst werden kann, habe ich die auto.arima()Funktion aus dem Prognosepaket verwendet, um herauszufinden, ob es ein AR (2) -Modell vorschlagen würde. Es stellt sich heraus, dass dieauto.arima() Funktion ein ARMA (2,2) -Modell vorschlägt; kein reines AR (2) Modell, aber das ist OK. Dies ist in Ordnung, da ein ARMA (2,2) -Modell eine AR (2) -Komponente enthält. Daher gilt dieselbe Regel (für stochastische Zyklen). Das heißt, wir können immer noch die oben genannte Bedingung überprüfen, um festzustellen, ob Sinuswellenvorhersagen erstellt werden.

Die Ergebnisse von auto.arima(y)sind unten gezeigt.

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

ϕ12+4ϕ2<01.73472+4(0.8324)<00.3202914<0
und wir stellen fest, dass die Bedingung tatsächlich erfüllt ist.

Das folgende Diagramm zeigt die Originalserie, y, die Passform des ARMA (2,2) -Modells und 14 Prognosen außerhalb der Stichprobe. Wie zu sehen ist, folgen die Vorhersagen außerhalb der Stichprobe einem Sinuswellenmuster.

Bildbeschreibung hier eingeben

Beachten Sie zwei Dinge. 1) Dies ist nur eine sehr schnelle Analyse (unter Verwendung eines automatisierten Tools), und eine ordnungsgemäße Behandlung würde die Befolgung der Box-Jenkins-Methodik erfordern. 2) ARIMA-Vorhersagen eignen sich gut für Kurzzeitvorhersagen. Sie können daher feststellen, dass Langzeitvorhersagen aus den Modellen in den Antworten von @David J. Harris und @Glen_b zuverlässiger sind.

Hoffentlich ist dies eine schöne Ergänzung zu einigen bereits sehr informativen Antworten.

Referenz : Vorhersage mit dynamischen Regressionsmodellen: Alan Pankratz, 1991 (John Wiley und Söhne, New York), ISBN 0-471-61528-5

Graeme Walsh
quelle
1

Die derzeitigen Methoden zum Anpassen einer Sinuskurve an einen bestimmten Datensatz erfordern ein erstes Erraten der Parameter, gefolgt von einem interaktiven Prozess. Dies ist ein nichtlineares Regressionsproblem. Eine andere Methode besteht darin, die nichtlineare Regression durch eine geeignete Integralgleichung in eine lineare Regression umzuwandeln. Dann ist keine erste Vermutung und kein iterativer Prozess erforderlich: Die Anpassung wird direkt erhalten. Für die Funktion y = a + r * sin (w * x + phi) oder y = a + b * sin (w * x) + c * cos (w * x) siehe Seiten 35-36 des Papiers "Régression sinusoidale", veröffentlicht auf Scribd: http://www.scribd.com/JJacquelin/documents Bei der Funktion y = a + p * x + r * sin (w * x + phi): Seite 49-51 des Kapitels "Gemischte lineare und sinusförmige Regressionen". Bei komplexeren Funktionen wird der allgemeine Vorgang im Kapitel "Generalisierte sinusförmige Regression" auf den Seiten 54-61 erläutert, gefolgt von einem Zahlenbeispiel y = r * sin (w * x + phi) + (b / x) + c * ln (x), Seiten 62-63

JJacquelin
quelle
0

Wenn Sie den niedrigsten und höchsten Punkt Ihrer Kosinus-Daten kennen, können Sie mit dieser einfachen Funktion alle Kosinus-Koeffizienten berechnen:

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)
}

Im Folgenden wird die Veränderung der Temperatur über den Tag mit einer Kosinusfunktion simuliert, indem die Stunden und Temperaturwerte für die niedrigste und wärmste Stunde eingegeben werden:

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")

Die Ausgabe ist unten: Kosinus berechnet aus niedrigsten und höchsten Punkten

IVIM
quelle
3
Dieser Ansatz scheint besonders empfindlich für zufällig aussehende Abweichungen von rein sinusförmigem Verhalten zu sein, so dass er auf fast alle Datensätze wie den in der Frage dargestellten nicht anwendbar wäre. Möglicherweise könnte es verwendet werden, um Startwerte für einige der anderen iterativen Ansätze bereitzustellen, die in diesem Thread vorgeschlagen werden.
Whuber
stimmen zu, es ist das einfachste, wäre gut für eine einfache Annäherung unter bestimmten Voraussetzungen
IVIM
0

Eine andere Option ist die Verwendung der generischen Funktion optim oder nls. Ich habe beide ausprobiert, keiner von ihnen ist völlig robust

Die folgenden Funktionen übernehmen die Daten in y und berechnen die 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])

Die Verwendung ist die folgende:

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)

Der folgende Code vergleicht die Daten

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)
NMech
quelle