Calculando a tendência de um ângulo quando ele cruza 360 -> 0

Aug 17 2020

Tenho uma variável que mede um ângulo que descreve a posição relativa de dois objetos (ou seja, pode variar de 0 a 359) e gostaria de quantificar como isso mudou ao longo do tempo.

Por exemplo, aqui temos a posição relativa dos dois itens variando 1 grau por ano:

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

No entanto, tomar a inclinação aqui não faz sentido por causa do "cruzamento" que ocorre no ano 2000. Tenho várias amostras diferentes, algumas com esse problema e outras não. Não sei a priori quais amostras terão esse problema, nem quando o crossover ocorrerá, então não posso simplesmente aplicar algum tipo de compensação (ou seja, adicionar 360 aos últimos 20 anos).

Existe uma maneira aceita de calcular tendências angulares, considerando o fato de que 0 = 360?

Respostas

1 whuber Aug 17 2020 at 23:00

Pense no ângulo$y$a qualquer momento$t$como o acúmulo de pequenas mudanças no ângulo. Simbolicamente, quando$f(t)$é a taxa de variação do ângulo no tempo$t$e$t_0$é o início das observações,

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

seu problema é esse$y(t)$foi gravado modulo$360$graus - talvez com algum erro$\epsilon(t).$Ou seja, você observou apenas os valores

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

Você pode, no entanto, reconstruir$y(t) + \epsilon(t)$desde que você tenha observações suficientemente frequentes. Por tempos sucessivos$t \lt s,$perceber

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

Onde$\delta$é igual à contribuição dos erros$\epsilon(s)-\epsilon(t)$ mais, talvez, algum múltiplo inteiro de$360$sempre que houver uma quebra angular entre$y^{*}(t)$e$y^{*}(s).$Agora, desde que o tamanho do erro total$|\epsilon(s)-\epsilon(t)|$é menos do que$180$graus e desde que o ângulo não gire mais de uma vez, podemos descobrir se ocorreu uma quebra: se$|\epsilon(s)-\epsilon(t)| \gt 180,$adicionar ou subtrair$360$graus de$\delta$colocá-lo no intervalo de$-180$para$+180$graus.

Embora não possamos observar esses erros diretamente, se estivermos amostrando com frequência suficiente para fazer os incrementos$y(t_i) - y(t_{i-1})$bastante pequeno, simplesmente aplicamos esse ajuste às diferenças observadas. Desta forma,

Em qualquer momento$|y^{*}(s)-y^{*}(t)| \gt 180,$adicionar ou subtrair$360$graus de$\delta$colocá-lo no intervalo de$-180$para$+180$graus.

Equivalentemente, calcule as diferenças módulo$180$mas expressá-los no intervalo de$-180$para$+180$graus em vez de (como é convencional) a faixa de$0$para$360.$

Vamos chamar o valor ajustado$\delta^{*}(t,s),$de modo a

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

Isso é igualdade, não igualdade módulo$360.$ Podemos agora remover o efeito de gravar os ângulos módulo$360$somando essas diferenças ajustadas. Quando as observações são feitas às vezes$t_0 \lt t_1\lt \cdots \lt t_n,$temos

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

O problema com módulo de computação$360$se foi: agora você pode usar qualquer procedimento que desejar para modelar a resposta$y^{*}(t).$


Aqui está uma ilustração com um conjunto de dados bastante difícil. Os dados foram gerados de acordo com o modelo$y(t) = 30t \mod 360$e observado anualmente de 1980 a 2020 com iid Erro normalmente distribuído do desvio padrão$60$graus (uma grande quantidade).

A tendência é quase imperceptível nos dados brutos, mas o algoritmo de ajuste de ângulo os alinhou visivelmente. Podemos ajustar um modelo de mínimos quadrados aos dados ajustados, por exemplo, produzindo este resultado:

A escala vertical expandida para os dados brutos mostra detalhes do ajuste e seus desvios. Aliás, neste exemplo, a estimativa da inclinação é$28.0 \pm 0.74$graus, não notavelmente diferente do verdadeiro valor de$30$graus (o valor-p para esta comparação é$1.1\%$).

Terminarei observando que quando o desvio padrão dos erros$\epsilon(t)$é grande (maior que$180/2/\sqrt{2} \approx 64$graus, aproximadamente), às vezes o ajuste angular ficará incorreto. Isso aparecerá nos resíduos do modelo como uma mudança repentina em um valor em torno de 360 ​​graus. Assim, uma análise rotineira dos resíduos do modelo pode detectar tais problemas, permitindo modificar os valores ajustados para um melhor ajuste. Os detalhes disso dependerão do seu modelo e do procedimento de adaptação.


Este Rcódigo criou as figuras. Em "ajustar os ângulos" mostra como o ajuste do ângulo pode ser calculado de forma eficiente.

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