Calcolo dell'andamento di un angolo quando attraversa 360 -> 0

Aug 17 2020

Ho una variabile che misura un angolo che descrive la posizione relativa di due oggetti (cioè può variare da 0 a 359) e vorrei quantificare come questo è cambiato nel tempo.

Ad esempio, qui abbiamo la posizione relativa dei due elementi che cambia di 1 grado all'anno:

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

Tuttavia, prendere la pendenza qui non ha senso a causa del "crossover" che si verifica nell'anno 2000. Ho un numero di campioni diversi, e alcuni hanno questo problema e altri no. Non so a priori quali campioni avranno questo problema, né quando si verificherà il crossover, quindi non posso semplicemente applicare una sorta di offset (ad esempio, aggiungere 360 ​​agli ultimi 20 anni).

Esiste un modo accettato per calcolare le tendenze angolari, tenendo conto del fatto che 0 = 360?

Risposte

1 whuber Aug 17 2020 at 23:00

Pensa all'angolo$y$in ogni momento$t$come l'accumulo di piccoli cambiamenti nell'angolo. Simbolicamente, quando$f(t)$è la velocità di variazione dell'angolo nel tempo$t$e$t_0$è l'inizio delle osservazioni,

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

Il tuo problema è quello$y(t)$è stato registrato modulo$360$gradi -- forse con qualche errore$\epsilon(t).$Cioè, hai osservato solo i valori

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

Puoi, tuttavia, ricostruire$y(t) + \epsilon(t)$a condizione di avere osservazioni sufficientemente frequenti. Per tempi successivi$t \lt s,$Avviso

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

dove$\delta$eguaglia il contributo degli errori$\epsilon(s)-\epsilon(t)$ più, forse, qualche multiplo intero di$360$ogni volta che c'è stata una rottura angolare tra$y^{*}(t)$e$y^{*}(s).$Ora, fornito la dimensione dell'errore totale$|\epsilon(s)-\epsilon(t)|$è meno di$180$gradi e ammesso che l'angolo non abbia girato più di una volta, possiamo capire se si è verificata un'interruzione: se$|\epsilon(s)-\epsilon(t)| \gt 180,$aggiungere o sottrarre$360$gradi da$\delta$per posizionarlo nell'intervallo da$-180$a$+180$gradi.

Sebbene non possiamo osservare direttamente questi errori, se stiamo campionando abbastanza frequentemente per fare gli incrementi$y(t_i) - y(t_{i-1})$abbastanza piccolo, applichiamo semplicemente questo aggiustamento alle differenze osservate. Così,

Ogni volta$|y^{*}(s)-y^{*}(t)| \gt 180,$aggiungere o sottrarre$360$gradi da$\delta$per posizionarlo nell'intervallo da$-180$a$+180$gradi.

In modo equivalente, calcola le differenze modulo$180$ma esprimili nell'intervallo da$-180$a$+180$gradi piuttosto che (come è convenzionale) l'intervallo da$0$a$360.$

Chiamiamo il valore regolato$\delta^{*}(t,s),$affinché

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

Questa è uguaglianza, non uguaglianza modulo$360.$ Possiamo ora rimuovere l'effetto della registrazione degli angoli modulo$360$sommando queste differenze aggiustate. Quando le osservazioni sono fatte a volte$t_0 \lt t_1\lt \cdots \lt t_n,$noi abbiamo

$$\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}$$

Il problema con il calcolo modulo$360$è andato: ora puoi usare qualsiasi procedura che ti piace per modellare la risposta$y^{*}(t).$


Ecco un'illustrazione con un set di dati abbastanza difficile. I dati sono stati generati secondo il modello$y(t) = 30t \mod 360$e osservato annualmente dal 1980 al 2020 con iid Errore normalmente distribuito della deviazione standard$60$gradi (una grande quantità).

La tendenza è appena percettibile nei dati grezzi, ma l'algoritmo di regolazione dell'angolo li ha allineati visibilmente. Possiamo adattare un modello dei minimi quadrati ai dati rettificati, ad esempio, producendo questo risultato:

La scala verticale espansa per i dati grezzi mostra i dettagli dell'adattamento e le loro deviazioni da esso. Per inciso, in questo esempio la stima della pendenza è$28.0 \pm 0.74$gradi, non notevolmente diverso dal vero valore di$30$gradi (il p-value per questo confronto è$1.1\%$).

Concluderò osservando che quando la deviazione standard degli errori$\epsilon(t)$è grande (maggiore di$180/2/\sqrt{2} \approx 64$gradi, approssimativamente), a volte la regolazione angolare sarà errata. Questo apparirà nei residui del modello come un cambiamento improvviso di un valore di circa 360 gradi. Pertanto, un'analisi di routine dei residui del modello può rilevare tali problemi, consentendo di modificare i valori regolati per un migliore adattamento. I dettagli di questo dipenderanno dal modello e dalla procedura di montaggio.


Questo Rcodice ha creato le figure. In "regola gli angoli" mostra come la regolazione dell'angolo può essere calcolata in modo efficiente.

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