- O sistema binário de raios-X GX 5-1
- Curvas de luz estelares da missão Kepler
adaptação do Time series analysis, Summer School in Statistics for Astronomers, Eric D. Feigelson, 2021
https://sites.psu.edu/astrostatistics/files/2021/02/Tutorial_4_Time_Series.html
Ver também Modern Statistical Methods for Astronomy with R Applications, E.D.Feigelson & G.J. Babu (2012)
GX 5-1 é um sistema binário que é forte emissor de raios-X. Nesse sistema, uma estrela normal transfere massa para uma estrela de neutrons via um disco de acresção, e os raios-X são emitidos na região mais interna, mais quente, desse disco.
O arquivo GX.dat contêm 65536 medidas de contagens de fotons em raios-X em bins igualmente espaçados por 1/128 segundos e foram obtidos pelo satélite japonês Ginga na década de 80.
As taxas de contagem são um processo poissoniano, mas como as contagens por bin são relativamente altas, podem ser aproximadas como um processo gaussiano.
leitura dos dados:
### Time series analysis
### Summer School in Statistics for Astronomers
### Eric D. Feigelson Spring 2021
# Read the dataset
GX.dat <- scan("GX.dat")
head(GX.dat)[1] 79 70 68 58 50 59summary(GX.dat) Min. 1st Qu. Median Mean 3rd Qu. Max.
30.00 62.00 69.00 68.82 76.00 116.00 Note que os dados serão colocados num objeto de classe ‘ts’ (“time series”), adequado para a análise de séries temporais em R:
# Convert to a class 'ts' (time series) R object
# tempo: note que a série temporal corresponde a 65536/128 = 512 segundos
GX.time <- seq(from=0, to=512, length.out=length(GX.dat))
GX.ts <- ts(GX.dat, GX.time) ; GX.ts.offset <- ts(GX.dat-30, GX.time)
# visualização dos dados:
plot.ts(GX.ts, ylab='GX 5-1 counts', xlab='Time (x 1/128 sec)', cex.lab=1.3, cex.axis=1.3, lwd=0.5,col='gold3')
Vamos comparar a distribuição de contagens com uma distribuição normal e com uma distribuição de Poisson.
# opções gráficas
options(jupyter.plot_scale=1)
options(repr.plot.width = 7, repr.plot.height = 5)
hist(GX.dat, breaks=100, xlim=c(40,100), ylim=c(0,3500), xlab='GX 5-1 counts',
font=2, font.lab=2, main='')
sdgauss = sd(GX.dat)
sdpoisson = sqrt(mean(GX.dat))
curve(dnorm(x,mean=mean(GX.dat), sd=sdgauss)*65536, lwd=3, add=T)
curve(dnorm(x,mean=mean(GX.dat), sd=sdpoisson)*65536, lwd=3, add=T, col='blue')
legend('topleft',legend=c('gaussiana','poisson'),col=c('black','blue'), text.col = c('black','blue'),bty='n')
O resultado não é consistente com um processo puramente poissoniano!
Vamos dar um zoom na série, vendo \(\sim10\)% dos dados:
plot.ts(GX.ts[1:6000], ylab='GX 5-1 counts', xlab='Time (x 1/128 sec)', cex.lab=1.3, cex.axis=1.3,col='cadetblue') 
será que tem algum comportamento semiperiódico presente?
Vamos exibir agora versões filtradas da série:
plot(GX.time,GX.dat, ylim=c(-10,115), xlab='Time (sec)', ylab='GX 5-1 counts',
cex.lab=1.3, cex.axis=1.3, type='n') # set up plot window but don't show any data
lines(ksmooth(GX.time, GX.dat+30, 'normal', bandwidth=7), lwd=2)
text(450, 110, 'Normal kernel') # Gaussian kernel density estimator with 7 bin FWHM bandwidth
lines(filter(GX.ts, sides=2, rep(1,7)/7), lwd=2)
text(450, 85, 'Moving average') # Moving average smoother with 7 bin bandwidth
lines(kernapply(GX.ts.offset, kernel('modified.daniell', 7)), lwd=2)
text(450, 50, 'Modified Daniell') # Moving average smoother with 1/2-weight at the end values of the span
lines(supsmu(GX.time, GX.dat-60, span=0.01), lwd=2)
text(400, 20, "Friedman's super-smoother") # A smoother with adaptive bandwidth from Friedman (1984)
lines(lowess(GX.time, GX.dat-80, 0.02), lwd=2)
text(400, 0, 'LOWESS local regression') # Cleveland's (1979) robust local polynomial smoother
Algumas das curvas acima sugerem variações com escalas de tempo de 20 a 50 segundos.
Uma ferramenta importante para analisar as séries temporais é a função de autocorrelação: ela analisa a correlação de uma variável entre instantes separados por um certo intervalo \(k\).
Dada uma série temporal \(\{x_1, x_2, ...\}\), a ACF(k) é a fração da variância total devido a valores com atraso de k intervalos: \[ ACF(k) = {\sum_{t=1}^{n-k} (x_t - \bar x)(x_{t-k}- \bar x) \over \sum_{t=1}^{n} (x_t - \bar x)^2} \]
par(mfrow=c(1,1))
acf(GX.dat)
# as linhas tracejadas azuis representam o intervalo de confiança de 95% A autocorrelação com k=1 é maior que 20%, o que mostra uma forte dependência do valor atual de valores anteriores.
Parece haver uma componente quasiperiódica com \(P \sim 6\) segundos (k=5).
A ACF mostra que essa estrutura é significativa, excedendo o intervalo de confiança de 95% (linhas tracejadas) até \(k \sim 30\).
Periodogramas são outra ferramenta importante na análise de séries temporais.
Periodograma de Schuster (clássico):
# Spectral periodogram
par(mfrow=c(1,1))
spec.pgram(GX.ts, log='no', main='')
Periodograma de Schuster suavisado:
# Spectral periodogram
par(mfrow=c(2,1))
spec.pgram(GX.ts, spans=50, log='no', main='', sub='')
spec.pgram(GX.ts, spans=200, taper=0.15, log='no', main='', sub='')
Os periodogramas mostram estruturas associadas a dois tipos de processos: o crescimento em direção a baixas frequências (abaixo de 0.05) é devido a ruído (red noise), e o pico alargado, em 0.17-0.23, é devido às oscilações quasiperiódicas.
Vamos agora considerar um modelo autoregressivo simples, onde o valor atual depende linearmente de um certo número de valores em tempos anteriores. Nesse caso, o algoritmo encontra que o melhor modelo (usando AIC) tem 27 coeficientes:
# Autoregressive modeling: AR
ARmod <- ar(GX.ts, method='ols')
print(ARmod)
Call:
ar(x = GX.ts, method = "ols")
Coefficients:
1 2 3 4 5 6 7 8
0.2047 0.0153 0.0046 0.0674 0.1086 0.0515 -0.0166 -0.0242
9 10 11 12 13 14 15 16
0.0059 0.0430 0.0370 -0.0053 -0.0097 -0.0073 0.0032 0.0241
17 18 19 20 21 22 23 24
-0.0010 -0.0064 -0.0167 0.0040 0.0177 0.0034 0.0025 -0.0067
25 26 27
-0.0050 0.0073 0.0113
Intercept: -4.655e-05 (0.03832)
Order selected 27 sigma^2 estimated as 96.18ARspec <- spec.ar(GX.ts, plot=F)
GXspec <- spec.pgram(GX.ts, span=101, main='', sub='', lwd=2)
lines(ARspec$freq, ARspec$spec, col='green', lwd=2)
text(0.4,450, cex=2, paste0('AR(', ARmod$order, ')'))
Uma forma de verificar se o modelo autoregressivo ajustado é bom, é calculando a função de autocorrelação dos resíduos do ajuste: espera-se que os resíduos não estejam correlacionados!
acf(na.omit(ARmod$resid), main='ACF dos resíduos')
Pode-se notar que o modelo ajusta muito bem a ACF, sem introduzir nenhuma componente periódica.
Existem vários pacotes que implementam modelos autoregressivos em R. Exemplos adicionais de aplicação a GX 5-1 podem ser encontrados em https://sites.psu.edu/astrostatistics/files/2021/02/Tutorial_4_Time_Series.html
A análise desta seção nos leva a concluir que GX 5-1 exibe variabilidade estocástica de baixa amplitude, com oscilações semiperiódicas de 20-50 segundos.
A missão Kepler mediu com alta precisão o brilho na faixa do visível de cerca de 150 mil estrelas cada 29.4 minutos por quase 4 anos, tendo detectado milhares de trânsitos exoplanetários, onde a passagem de um planeta em frente a uma estrela provoca uma diminuição de seu brilho.
Nessa seção vamos considerar um caso de variabilidade estelar que não tem nada a ver com exoplanetas e que ilustra as dificuldades que podem aparecer na busca por esses objetos.
Esse exemplo considera séries temporais igualmente espaçadas. Mas as séries observadas têm, em geral, muitos dados faltantes e os métodos que vamos considerar exigem séries completas. A solução é imputar dados, isto é, atribuir valores aos dados faltantes (NA, Not Available) usando algum procedimento estatístico.
Essa estrela mostra variabilidade quase-periódica devido à presença de rotação e manchas estelares.
# leitura dos dados: vetor pré-processado com medidas de fluxo igualmente espaçadas no tempo
Kepler2 <- read.table('Kepler2.dat')[[1]]
length(Kepler2)[1] 71427# fração de NAs:
length(which(is.na(Kepler2))) / length(Kepler2) [1] 0.2883783# ~29% NAs
# visualização
par(mfrow=c(1,1))
plot(Kepler2, type='l', xlab='Time')
Note que 29% dos dados são NAs. Vamos substituir os NAs por valores “imputados”.
Há inúmeras técnicas de imputação; vamos usar um filtro de Kalman:
Em estatística, o filtro de Kalman é um método matemático criado por Rudolf Kalman. Seu propósito é utilizar medições de grandezas realizadas ao longo do tempo (contaminadas com ruído e outras incertezas) e gerar resultados que tendam a se aproximar dos valores reais das grandezas medidas e valores associados. O filtro de Kalman apresenta diversas aplicações e é uma parte essencial do desenvolvimento de tecnologias espaciais e militares. (Wikipedia)
library(imputeTS)
Kepler2_impute <- na_kalman(Kepler2) # impute NAs with Kalman smootherWarning in stats::StructTS(data, ...): possible convergence problem: 'optim'
gave code = 52 and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'# veja o antes e o depois da imputação:
par(mfrow=c(2,1)) ; par(mar=c(5,4,1,2))
plot(Kepler2, type='l', xlim=c(60000,70000), xlab='Time', ylab='Kepler 2')
plot(Kepler2_impute, type='l', xlim=c(60000,70000), xlab='Time',
ylab='Kepler 2 imputed')
Vamos examinar esta curva de luz (imputada) com 2 periodogramas:
- análise de Fourier usando a função spec.pgram
- o periodograma de Lomb-Scargle, muito usado em astronomia
# periodograma com análise de Fourier
par(mfrow=c(1,1))
spec.pgram(Kepler2_impute, xlim=c(0,0.005), spans=5, taper=0.0,
main='', ylab='Fourier', sub='')
# periodograma de Lomb-Scargle
library(lomb)
lsp(Kepler2, from=0.00, to=0.005, ylab='Lomb-Scargle', main='')
Ambas as técnicas sugerem dois períodos principais na curva de luz deste objeto.
Há muitas outras técnicas que todos os interessados em trabalhar com séries temporais precisam conhecer!
A série do número de manchas dolares com o tempo é bastante complexa, como pode ser vislumbrado por esse exercício.
A tabela SN_m_tot_V2.0.txt contém o número mensal de manchas solares entre 1749 e maio de 2024 (download em 3/06/2024), como disponível em https://www.sidc.be/SILSO/infosnmtot
A tabela contém: o ano, o mês, o ano com a fração do ano no meio do mês correspondente, o número de manchas solares, o desvio padrão do número de manchas se disponível, o número de observações se disponível.
- Leia os dados da tabela e faça um gráfico do número de manchas solares em função do tempo
- Suavize essa figura usando um kernel gaussiano com largura de banda de 5 meses. Comente o resultado: como se compara com a curva do item 1?
- Faça um plot da função de autocorrelação. Como você interpreta o resultado? Atenção: para visualizar a acf considere lags de até 1000 meses (faça ?acf e use o parâmetro lag.max)
- Faça um periodograma de Lomb-Scargle. Comente sobre as frequências relevantes que você vê.
- Faça um modelo autoregressivo AR. Compare com o ajuste obtido acima para GX 5-1.