Menghitung tren sudut saat melintasi 360 -> 0

Aug 17 2020

Saya memiliki variabel yang mengukur sudut yang menggambarkan posisi relatif dari dua objek (yaitu, dapat berkisar dari 0-359), dan saya ingin mengukur bagaimana hal ini telah berubah seiring waktu.

Misalnya, di sini kita memiliki posisi relatif dari dua item yang berubah sebesar 1 derajat per tahun:

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

Namun, mengambil kemiringan di sini tidak ada artinya karena "persilangan" yang terjadi pada tahun 2000. Saya memiliki sejumlah sampel yang berbeda, dan beberapa memiliki masalah ini dan beberapa tidak. Saya tidak tahu apriori sampel mana yang akan mengalami masalah ini, atau kapan persilangan terjadi, jadi saya tidak bisa hanya menerapkan semacam offset (yaitu, menambahkan 360 ke 20 tahun terakhir).

Apakah ada cara yang diterima untuk menghitung tren sudut, dengan memperhitungkan fakta bahwa 0 = 360?

Jawaban

1 whuber Aug 17 2020 at 23:00

Pikirkan sudutnya $y$ kapan saja $t$sebagai akumulasi perubahan kecil di sudut. Secara simbolis, kapan$f(t)$ adalah laju perubahan sudut pada saat itu $t$ dan $t_0$ adalah awal dari pengamatan,

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

Masalahmu adalah itu $y(t)$ telah direkam modulo $360$ derajat - mungkin dengan beberapa kesalahan $\epsilon(t).$ Artinya, Anda hanya mengamati nilainya

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

Anda bisa, bagaimanapun, merekonstruksi $y(t) + \epsilon(t)$asalkan Anda memiliki pengamatan yang cukup sering. Untuk waktu yang berurutan$t \lt s,$ memperhatikan

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

dimana $\delta$ sama dengan kontribusi kesalahan $\epsilon(s)-\epsilon(t)$ ditambah, mungkin, beberapa kelipatan integral dari$360$ setiap kali ada pemutusan sudut antara $y^{*}(t)$ dan $y^{*}(s).$Sekarang, asalkan ukuran total error$|\epsilon(s)-\epsilon(t)|$ kurang dari $180$derajat dan asalkan sudutnya tidak berputar lebih dari sekali, kita dapat mengetahui apakah terjadi penembusan: jika$|\epsilon(s)-\epsilon(t)| \gt 180,$ menambah atau mengurangi $360$ derajat dari $\delta$ untuk menempatkannya ke dalam jeda dari $-180$ untuk $+180$ derajat.

Meskipun kami tidak dapat mengamati kesalahan ini secara langsung, jika kami mengambil sampel cukup sering untuk melakukan peningkatan $y(t_i) - y(t_{i-1})$cukup kecil, kami hanya menerapkan penyesuaian ini pada perbedaan yang diamati. Jadi,

Kapanpun $|y^{*}(s)-y^{*}(t)| \gt 180,$ menambah atau mengurangi $360$ derajat dari $\delta$ untuk menempatkannya ke dalam jeda dari $-180$ untuk $+180$ derajat.

Secara ekivalen, hitunglah perbedaan modulo $180$ tetapi ekspresikan dalam rentang dari $-180$ untuk $+180$ derajat daripada (seperti konvensional) kisaran dari $0$ untuk $360.$

Mari kita sebut nilai yang disesuaikan $\delta^{*}(t,s),$ yang seperti itu

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

Ini adalah kesetaraan, bukan kesetaraan modulo$360.$ Kami sekarang dapat menghapus efek perekaman sudut modulo $360$dengan menjumlahkan perbedaan yang disesuaikan ini. Saat observasi dilakukan pada waktu-waktu tertentu$t_0 \lt t_1\lt \cdots \lt t_n,$ kita punya

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

Masalah dengan modulo komputasi $360$hilang: Anda sekarang dapat menggunakan prosedur apa pun yang Anda suka untuk memodelkan respons$y^{*}(t).$


Berikut ilustrasi dengan dataset yang cukup sulit. Data dihasilkan sesuai model$y(t) = 30t \mod 360$dan diamati setiap tahun dari tahun 1980 hingga 2020 dengan iid Biasanya didistribusikan kesalahan standar deviasi$60$ derajat (jumlah besar).

Trennya hampir tidak terlihat dalam data mentah, tetapi algoritme penyesuaian sudut tampak menyelaraskannya. Kami dapat memasukkan model kuadrat terkecil ke data yang disesuaikan, misalnya, menghasilkan hasil ini:

Skala vertikal yang diperluas untuk data mentah menunjukkan detail kecocokan dan penyimpangannya. Secara kebetulan, dalam contoh ini perkiraan kemiringannya adalah$28.0 \pm 0.74$ derajat, tidak terlalu berbeda dari nilai sebenarnya $30$ derajat (nilai p untuk perbandingan ini adalah $1.1\%$).

Saya akan mengakhiri dengan mengatakan bahwa ketika standar deviasi kesalahan $\epsilon(t)$ besar (lebih besar dari $180/2/\sqrt{2} \approx 64$derajat, secara kasar), terkadang penyesuaian sudut akan salah. Ini akan muncul dalam residual model sebagai perubahan mendadak dengan nilai sekitar 360 derajat. Dengan demikian, analisis rutin dari residual model dapat mendeteksi masalah tersebut, memungkinkan Anda untuk memodifikasi nilai yang disesuaikan agar lebih sesuai. Detailnya tergantung pada model dan prosedur pemasangan Anda.


RKode ini menciptakan angka-angka. Pada "sesuaikan sudut" ini menunjukkan bagaimana penyesuaian sudut dapat dihitung secara efisien.

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