Ajuster un terme sinusoïdal aux données

26

Bien que j'ai lu ce post, je n'ai toujours aucune idée de comment l'appliquer à mes propres données et j'espère que quelqu'un pourra m'aider.

J'ai les données suivantes:

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

Et maintenant je veux simplement rentrer dans une onde sinusoïdale

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

avec les quatre inconnues , , et à elle.AωϕC

Le reste de mon code est le suivant

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

Mais le résultat est vraiment médiocre.

Ajustement sinusoïdal

J'apprécierais beaucoup toute aide.

À votre santé.

Pascal
la source
Vous essayez d'adapter une onde sinusoïdale aux données ou essayez-vous d'adapter une sorte de modèle harmonique avec un composant sinusoïdal et cosinusoïdal? Il y a une fonction harmonique dans le package TSA dans R que vous voudrez peut-être vérifier. Ajustez votre modèle en utilisant cela et voyez quel type de résultats vous obtenez.
Eric Peterson
5
Avez-vous essayé différentes valeurs de départ? Votre fonction de perte n'est pas convexe, donc différentes valeurs de départ peuvent conduire à des solutions différentes.
Stefan Wager
1
Dites-nous en plus sur les données. Il existe généralement une périodicité connue, de sorte qu'il n'est pas nécessaire d'estimer les données. Est-ce une série chronologique ou autre chose? C'est beaucoup plus facile si vous pouvez ajuster des termes sinus et cosinus séparés par un modèle linéaire.
Nick Cox
2
Le fait d'avoir une période inconnue rend votre modèle non linéaire (un tel événement est mentionné dans la réponse sélectionnée au message lié). Étant donné que, les autres paramètres sont conditionnellement linéaires; pour certaines routines LS non linéaires, ces informations sont importantes et peuvent améliorer le comportement. Une option pourrait être d'utiliser des méthodes spectrales pour obtenir la période et la condition à ce sujet; une autre serait de mettre à jour la période et les autres paramètres via une optimisation non linéaire et linéaire respectivement de manière itérative.
Glen_b -Reinstate Monica
(Je viens de modifier la réponse ici pour faire du cas particulier d'une période inconnue un exemple explicite de ce qui peut le rendre non linéaire.)
Glen_b -Reinstate Monica

Réponses:

18

Si vous voulez juste une bonne estimation de et que vous ne vous souciez pas beaucoup de son erreur standard:ω

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

tracé sinus

(Un meilleur ajustement expliquerait peut-être encore les valeurs aberrantes de cette série, en réduisant leur influence.)

---

Si vous voulez avoir une idée de l'incertitude dans , vous pouvez utiliser la vraisemblance du profil ( pdf1 , pdf2 - les références pour obtenir des CI ou SE approximatifs à partir de la vraisemblance du profil ou de ses variantes ne sont pas difficiles à localiser)ω

(Alternativement, vous pouvez alimenter ces estimations dans nls ... et les démarrer déjà convergées.)

Glen_b -Reinstate Monica
la source
(+1) belle réponse. J'ai essayé d'adapter le modèle linéaire avec lm(y~sin(2*pi*t)+cos(2*pi*t)mais cela n'a pas fonctionné (le costerme était toujours 1). Juste par curiosité: que font les deux premières raies (je sais qui spectrumestime la densité spectrale)?
COOLSerdash
1
t2*pi*t
1
@COOLSerdash (ctd) - La 2e ligne trouve la fréquence associée au plus grand pic du spectre et inverse pour identifier la période. Au moins dans ce cas (mais je soupçonne plus largement), les valeurs par défaut identifient essentiellement la période qui maximise la probabilité si étroitement que j'ai supprimé les étapes que j'ai eues pour maximiser la probabilité de profil dans la région autour de cette période. La fonction specdans TSA peut être meilleure (elle semble avoir plus d'options, dont l'une peut parfois être importante), mais dans ce cas, le pic principal était exactement au même endroit que pour, spectrumdonc je n'ai pas pris la peine.
Glen_b -Reinstate Monica
@Glen_b cette méthode fait des merveilles pour mon cas d'utilisation. J'ai également besoin d'ajuster une courbe cos (x), mais cela ne fonctionne pas aussi bien ... J'ai changé le reslmen reslm <- lm(y ~ cos(2*pi/per*t)+tan(2*pi/per*t))mais cela ne semble pas correct. des indices?
Amit Kohli
Pourquoi avez-vous un terme de bronzage là-dedans?
Glen_b -Reinstate Monica
15

2π/20

Quand j'ai mis cela dans nlsla startliste de, j'ai obtenu une courbe qui était beaucoup plus raisonnable, bien qu'elle ait encore des biais systématiques.

En fonction de votre objectif avec cet ensemble de données, vous pouvez essayer d'améliorer l'ajustement en ajoutant des termes supplémentaires ou en utilisant une approche non paramétrique comme un processus gaussien avec un noyau périodique.

Ajustement sinusoïdal

Choisir automatiquement une valeur de départ

Si vous souhaitez choisir la fréquence dominante, vous pouvez utiliser une transformée de Fourier rapide (FFT). C'est loin de mon domaine d'expertise, donc je vais laisser d'autres personnes remplir les détails s'ils le souhaitent (en particulier sur les étapes 2 et 3), mais le Rcode ci-dessous devrait fonctionner.

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

Vous pouvez également tracer abs(truncated.fft)pour voir s'il existe d'autres fréquences importantes, mais vous devrez jouer un peu avec la mise à l'échelle de l'axe des x.

De plus, je crois que @Glen_b a raison que le problème est convexe une fois que vous connaissez les oméga (ou peut-être avez-vous aussi besoin de connaître le phi? Je ne suis pas sûr). Dans tous les cas, connaître les valeurs de départ pour les autres paramètres ne devrait pas être aussi important que pour les oméga s'ils sont dans le bon stade. Vous pourriez probablement obtenir des estimations décentes des autres paramètres de la FFT, mais je ne sais pas comment cela fonctionnerait.

David J. Harris
la source
1
Merci pour cet indice. Juste pour clarifier un peu: les données font partie d'un microréseau dans lequel la périodicité des gènes a été mesurée dans le temps, c'est-à-dire que les données présentées sont les données d'expression d'un gène. Le problème est maintenant que je veux appliquer cette méthode à environ 40k gènes ayant tous des périodicités et des amplitudes différentes. Il est donc très important de trouver un bon ajustement indépendamment des conditions initiales.
Pascal
1
@Pascal Voir mes mises à jour ci-dessus pour une recommandation pour choisir automatiquement la valeur de départ pour les oméga.
David J. Harris
2
ϕuneb
Je me demande où les valeurs x entrent en jeu ici. Bien sûr, cela fait une différence pour les oméga, que les valeurs y données soient séparées par 1 ou par 5 étapes x, n'est-ce pas?
knub
1
Astuce de programmation non liée à la question: prudence lors de la désignation des objets R comme foo.bar. Cela est dû à la façon dont R spécifie les méthodes pour les classes .
Firebug
10

Comme alternative à ce qui a déjà été dit, il peut être intéressant de noter qu'un modèle AR (2) de la classe des modèles ARIMA peut être utilisé pour générer des prévisions avec un modèle d'onde sinusoïdale.

yt=C+ϕ1yt-1+ϕ2yt-2+unet
Cϕ1ϕ2unet

ϕ12+4ϕ2<0.

Panratz (1991) nous apprend ce qui suit sur les cycles stochastiques:

Un modèle de cycle stochastique peut être considéré comme un modèle d'onde sinusoïdale déformé dans le modèle de prévision: il s'agit d'une onde sinusoïdale avec une période, une amplitude et un angle de phase stochastiques (probabilistes).

Pour voir si un tel modèle pouvait être ajusté aux données, j'ai utilisé la auto.arima()fonction du package de prévision pour savoir s'il suggérerait un modèle AR (2). Il s'avère que la auto.arima()fonction suggère un modèle ARMA (2,2); pas un pur modèle AR (2), mais c'est OK. C'est OK car un modèle ARMA (2,2) contient un composant AR (2), donc la même règle (sur les cycles stochastiques) s'applique. Autrement dit, nous pouvons toujours vérifier la condition susmentionnée pour voir si des prévisions d'ondes sinusoïdales seront produites.

Les résultats de auto.arima(y)sont présentés ci-dessous.

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)<0-0,3202914<0

Le graphique ci-dessous montre la série d'origine, y, l'ajustement du modèle ARMA (2,2) et 14 prévisions hors échantillon. Comme on peut le voir, les prévisions hors échantillon suivent un schéma sinusoïdal.

entrez la description de l'image ici

Gardez à l'esprit deux choses. 1) Il s'agit simplement d'une analyse très rapide (à l'aide d'un outil automatisé) et un traitement approprié impliquerait de suivre la méthodologie de Box-Jenkins. 2) Les prévisions ARIMA sont bonnes pour les prévisions à court terme, vous pouvez donc trouver que les prévisions à long terme des modèles dans les réponses de @David J. Harris et @Glen_b sont plus fiables.

Enfin, j'espère que c'est un bon ajout à certaines réponses déjà très instructives.

Référence : Prévision avec des modèles de régression dynamique: Alan Pankratz, 1991, (John Wiley and Sons, New York), ISBN 0-471-61528-5

Graeme Walsh
la source
1

Les méthodes actuelles pour adapter une courbe sin à un ensemble de données donné nécessitent une première estimation des paramètres, suivie d'un processus interactif. Il s'agit d'un problème de régression non linéaire. Une autre méthode consiste à transformer la régression non linéaire en une régression linéaire grâce à une équation intégrale pratique. Ensuite, il n'y a pas besoin de conjecture initiale et pas besoin de processus itératif: l'ajustement est directement obtenu. Dans le cas de la fonction y = a + r * sin (w * x + phi) ou y = a + b * sin (w * x) + c * cos (w * x), voir pages 35-36 du document "Régression sinusoidale" publiée sur Scribd: http://www.scribd.com/JJacquelin/documents Dans le cas de la fonction y = a + p * x + r * sin (w * x + phi): pages 49-51 du chapitre "Régressions mixtes linéaires et sinusoïdales". Dans le cas de fonctions plus compliquées, le processus général est expliqué dans le chapitre "Régression sinusoïdale généralisée" pages 54-61, suivi d'un exemple numérique y = r * sin (w * x + phi) + (b / x) + c * ln (x), pages 62-63

JJacquelin
la source
0

Si vous connaissez le point le plus bas et le plus haut de vos données d'aspect cosinus, vous pouvez utiliser cette fonction simple pour calculer tous les coefficients cosinus:

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

Ci-dessous, il est utilisé pour simuler la variation de la température tout au long de la journée avec une fonction cosinus, en entrant les heures et les valeurs de température pour l'heure la plus basse et la plus chaude:

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

La sortie est ci-dessous: Cosinus calculé à partir des points les plus bas et les plus hauts

IVIM
la source
3
Cette approche semblerait être particulièrement sensible à toute dérogation aléatoire au comportement sinusoïdal pur, ce qui la rendrait inapplicable à presque tous les ensembles de données comme celui illustré dans la question. En théorie, il pourrait être utilisé pour fournir des valeurs de départ pour certaines des autres approches itératives suggérées dans ce fil.
whuber
d'accord, c'est le plus simple, serait bon pour une simple approximation sous certaines hypothèses
IVIM
0

Une autre option utilise la fonction générique optim ou nls. J'ai essayé les deux aucun d'eux n'est complètement robuste

Les fonctions suivantes prennent les données en y et calculent les paramètres.

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

l'utilisation est la suivante:

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)

Le code suivant compare les données

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
la source