Calcul de la tendance d'un angle lorsqu'il passe par 360 -> 0

Aug 17 2020

J'ai une variable qui mesure un angle décrivant la position relative de deux objets (c'est-à-dire, peut varier de 0 à 359), et je voudrais quantifier comment cela a changé au fil du temps.

Par exemple, nous avons ici la position relative des deux éléments changeant de 1 degré par an :

year <- seq(1981, 2020)
angle <- c(seq(341, 359), seq(0, 20))

Cependant, prendre la pente ici n'a aucun sens à cause du "croisement" qui se produit en l'an 2000. J'ai un certain nombre d'échantillons différents, et certains ont ce problème et d'autres non. Je ne sais pas a priori quels échantillons auront ce problème, ni quand le croisement se produit, donc je ne peux pas simplement appliquer une sorte de décalage (c'est-à-dire ajouter 360 aux 20 dernières années).

Existe-t-il une manière acceptée de calculer les tendances angulaires, en tenant compte du fait que 0 = 360 ?

Réponses

1 whuber Aug 17 2020 at 23:00

Pensez à l'angle$y$à tout moment$t$que l'accumulation de petits changements dans l'angle. Symboliquement, quand$f(t)$est le taux de variation de l'angle au moment$t$et$t_0$est le début des observations,

$$y(t) = y(t_0) + \int_{t_0}^t f(t)\,\mathrm{d}t.$$

Votre problème est que$y(t)$a été enregistré modulo$360$degrés - peut-être avec une erreur$\epsilon(t).$Autrement dit, vous n'avez observé que les valeurs

$$y^{*}(t) = y(t) + \epsilon(t) \mod 360.$$

Vous pouvez cependant reconstruire$y(t) + \epsilon(t)$à condition d'avoir des observations suffisamment fréquentes. Pour des temps successifs$t \lt s,$remarquer

$$y^{*}(s) - y^{*}(t) = y(s) - y(t) + \epsilon(s) - \epsilon(t) \mod 360 = \int_t^s f(t)\,\mathrm{d}t + \delta$$

$\delta$est égal à la contribution des erreurs$\epsilon(s)-\epsilon(t)$ plus, peut-être, un multiple entier de$360$chaque fois qu'il y a eu une rupture angulaire entre$y^{*}(t)$et$y^{*}(s).$Maintenant, à condition que la taille de l'erreur totale$|\epsilon(s)-\epsilon(t)|$est inférieur à$180$degrés et à condition que l'angle n'ait pas tourné plus d'une fois, nous pouvons déterminer si une rupture s'est produite : si$|\epsilon(s)-\epsilon(t)| \gt 180,$ajouter ou soustraire$360$degrés de$\delta$pour le placer dans l'intervalle de$-180$à$+180$degrés.

Bien que nous ne puissions pas observer ces erreurs directement, si nous échantillonnons assez fréquemment pour effectuer les incréments$y(t_i) - y(t_{i-1})$assez faible, nous appliquons simplement cet ajustement aux différences observées. Ainsi,

À tout moment$|y^{*}(s)-y^{*}(t)| \gt 180,$ajouter ou soustraire$360$degrés de$\delta$pour le placer dans l'intervalle de$-180$à$+180$degrés.

De manière équivalente, calculez les différences modulo$180$mais exprimez-les dans la gamme de$-180$à$+180$degrés plutôt que (comme c'est conventionnel) la gamme de$0$à$360.$

Appelons la valeur ajustée$\delta^{*}(t,s),$pour que

$$y^{*}(s) - y^{*}(t) = \int_t^s f(t)\,\mathrm{d}t + \delta(t,s)^{*}.$$

C'est l' égalité, pas l'égalité modulo$360.$ Nous pouvons maintenant supprimer l'effet de l'enregistrement des angles modulo$360$en additionnant ces différences ajustées. Lorsque les observations sont faites à des moments$t_0 \lt t_1\lt \cdots \lt t_n,$Nous avons

$$\begin{aligned} y^{*}(t_i) &= y^{*}(t_0) + \left[y^{*}(t_1) - y^{*}(t_0)\right] + \cdots + \left[y^{*}(t_i) - y^{*}(t_{i-1})\right] \\ &=y(t_0) + \int_{t_0}^{t_i} f(t)\,\mathrm{d}t + \delta(t_0,t_1)^{*} + \delta(t_1,t_2)^{*} + \cdots + \delta(t_{i-1},t_i)^{*} \\ &= y(t_i) + \left[\epsilon(t_i) - \epsilon(t_0)\right]. \end{aligned}$$

Le problème du calcul modulo$360$est parti : vous pouvez maintenant utiliser la procédure de votre choix pour modéliser la réponse$y^{*}(t).$


Voici une illustration avec un jeu de données assez difficile. Les données ont été générées selon le modèle$y(t) = 30t \mod 360$et observé annuellement de 1980 à 2020 avec iid Erreur normalement distribuée de l'écart type$60$degrés (une grande quantité).

La tendance est à peine perceptible dans les données brutes, mais l'algorithme d'ajustement de l'angle les a visiblement alignées. Nous pouvons ajuster un modèle des moindres carrés aux données ajustées, par exemple, en produisant ce résultat :

L'échelle verticale développée pour les données brutes montre les détails de l'ajustement et leurs écarts par rapport à celui-ci. Incidemment, dans cet exemple, l'estimation de la pente est$28.0 \pm 0.74$degrés, pas remarquablement différent de la vraie valeur de$30$degrés (la valeur p pour cette comparaison est$1.1\%$).

Je terminerai en remarquant que lorsque l'écart type des erreurs$\epsilon(t)$est grand (supérieur à$180/2/\sqrt{2} \approx 64$degrés, à peu près), parfois le réglage angulaire sera incorrect. Cela apparaîtra dans les résidus du modèle comme un changement soudain d'une valeur d'environ 360 degrés. Ainsi, une analyse de routine des résidus du modèle peut détecter de tels problèmes, vous permettant de modifier les valeurs ajustées pour un meilleur ajustement. Les détails de cela dépendront de votre modèle et de la procédure de montage.


Ce Rcode a créé les chiffres. À "ajuster les angles", il montre comment le réglage de l'angle peut être calculé efficacement.

#
# Specify the data-generation process.
#
year <- 1980:2020 # Dates to use
beta <- 30        # Annual rate of change
sigma <- 60       # Error S.D.
#
# Generate the data.
#
set.seed(17)
angle <- (year * beta + rnorm(length(year), 0, sigma)) %% 360
X <- data.frame(year, angle)
#
# Adjust the angles.
#
X$`total angle` <- with(X, {
  d <- (diff(angle) + 180) %% 360 - 180
  cumsum(c(angle[1], d))
})
#
# Fit a model to the adjusted angles.
#
fit <- lm(`total angle` ~ year, X)
#
# Analyze the fit.
#
b <- coefficients(fit)
y.hat <- predict(fit)

#--Compute dates the fit must wrap around from 360 to 0:
y.breaks <- seq(floor(min(y.hat) / 360)*360, max(y.hat), by=360)
year.breaks <- (y.breaks - b[1]) / b[2]

#--Make the plots:
u <- ceiling(max(X$`total angle`)/360)
par(mfcol=c(1,2))

#--The fits:
plot(X$year, X$angle, pch=19, ylim=c(0, 360), yaxp=c(0, 360, 4),
     col="gray", ylab="Angle (degrees)", xlab="Year",
     main="Raw Data and Fit")
for (x in year.breaks) 
  abline(c(-x * b[2], b[2]), col="Red", lwd=2)

plot(X$year, X$`total angle`, ylim=c(0,u*360),  yaxp=c(0, u*360, u),
     xlab="Year", ylab="Total angle",
     main="Adjusted Data and Fit")
abline(fit, col="Red", lwd=2)

#--The raw data:
plot(X$year, X$angle, ylim=c(0,u*360),  yaxp=c(0, u*360, u),
     pch=19, col="gray", ylab="Angle (degrees)", xlab="Year",
     main="Raw Data")

plot(X$year, X$`total angle`, ylim=c(0,u*360),
     yaxp=c(0, u*360, u),
     xlab="Year", ylab="Total angle",
     main="Adjusted Data")
par(mfcol=c(1,1))