tópicos:

  1. bootstrap
  1. testes de hipótese
  • 2.1 comparação de médias e variâncias: a distância à LMC
  • 2.2 testes com simulações de bootstrap
  • 2.3 Einstein x Newton
  • 2.4 teste de proporção
  • 2.5 exercícios com a biblioteca MASS do R
  1. comparação de distribuições
  • 3.1 comparação de distribuições de variáveis nominais com o \(\chi^2\)
  • 3.2 o teste de Kolmogorov- Smirnov (KS)
  1. análise de correlação
  1. o caso da galáxia sem matéria escura

exercícios

1. bootstrap

O bootstrap é um procedimento de simulação baseado nos dados disponíveis.

Vamos ilustrar o bootstrap considerando o ajuste de uma reta: \[ y = a+bx\] onde \(w=(a,b)\) são os parâmetros do modelo

Vamos usar o bootstrap para estimar os erros dos parâmetros

Vamos inicialmente gerar um modelo de reta e adicionar erros gaussianos

# semente para assegurar a reprodutibilidade dos resultados
set.seed(12345)

# simulação dos dados
# n valores de x: entre 0 e 15
n = 10
x=runif(n,min=0,max=15)

# parâmetros da reta
afit=1
bfit=2

# variável dependente
y=afit+bfit*x

# erros: desvio padrão 1.8
err=rnorm(n,mean=0,sd=1.8)
y=y+err

# visualização dos dados
par(mfrow = c(1,1))
plot(x,y,main='y=1+2x     sd=1.8',pch=20, col='red',cex=2)
abline(1,2,lwd=2)

Vamos ajustar uma reta aos dados usando a função lm():

# ajuste de uma reta com a função lm() (*linear model*)
ajuste = lm(y~x)
summary(ajuste)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.96520 -0.88621  0.00849  1.08039  2.78521 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.137      1.563   1.367    0.209    
## x              1.866      0.151  12.360 1.71e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.814 on 8 degrees of freedom
## Multiple R-squared:  0.9502, Adjusted R-squared:  0.944 
## F-statistic: 152.8 on 1 and 8 DF,  p-value: 1.711e-06
# "ajuste" é uma variável tipo lista, cujos valores podem ser estimados via ajuste$variável
# os valores dos parâmetros são
# intercepto
ajuste$coefficients[1]
## (Intercept) 
##    2.137335
# inclinação
ajuste$coefficients[2]
##        x 
## 1.866421
# note que a função lm() também dá os erros nos parâmetros
# neste exemplo o erro no intercepto é muito maior que o na inclinação 
# veja também o valor p associado à estimativa de cada parâmetro: enquanto o valor p do intercepto é muito alto, o da inclinação é muito menor

# visualização do resultado
par(mfrow = c(1,1))
# plot dos pontos
plot(x,y,main='modelo:  y=1+2x,  sd=1.8',pch=20, col='red',cex=2)
# linha teórica
abline(1,2,lwd=3)  #a=1, b=2
# linha ajustada
abline(ajuste$coefficients[1],ajuste$coefficients[2],lwd=3,col='yellow3')
legend('topleft',,legend=c('modelo','ajuste'),col=c('black','yellow3'),  text.col = c('black','yellow3'),bty="n")

Vamos estimar os erros nos parâmetros fazendo agora um grande número de simulações de bootstrap; em cada simulação:

  • geramos um novo conjunto de dados por amostragem com substituição
  • ajustamos uma reta a esses dados e guardamos o valor dos parâmetros da reta
# para medir o tempo de processamento desse módulo em meu laptop
t0 = Sys.time()

# distribuição de a e b por bootstrap
# número de simulações
nboot=10000

# crio uma matriz para armazenar os 2 parâmetros de cada uma das nboot simulações
w=matrix(0,ncol=2,nrow=nboot)
# d0 é o vetor com os n=10 índices dos dados que serão amostrados
d0=1:n

# simulações
for (i in 1:nboot) {
# resampling
# sample() cria um novo vetor de dados amostrando d0 com substituição
dboot=sample(d0,n, replace = TRUE)
xl=x[dboot]
yl=y[dboot]
#ajuste = lm(yl~xl,weights = 1/err^2)
ajuste = lm(yl~xl)
w[i,]=ajuste$coefficients
}

# tempo de processamento
Sys.time()-t0
## Time difference of 5.058596 secs

Agora podemos analisar os resultados:

# parâmetros das simulações
asim=w[,1]
bsim=w[,2]

# média dos parâmetros
mean(asim)
## [1] 2.117057
mean(bsim)
## [1] 1.865941
# desvio padrão dos parâmetros:
sd(asim)
## [1] 1.817053
sd(bsim)
## [1] 0.158168
# visualização
par(mfrow = c(1,2))
hist(asim,xlab='a',main='10000 simulações',breaks=50,col='orange')
abline(v=mean(asim),lwd=2)
hist(bsim,xlab='b',main='y=a+bx',breaks=50,col='orange')
abline(v=mean(bsim),lwd=2)

as linhas pretas são os valores médios das simulações

compare os erros nos parâmetros (desvios padrão) calculados por bootstrap com os calculados por lm()

em geral se usa bootstrap para estimar erros nos parâmetros; as estimativas dos parâmetros vêm do ajuste (por exemplo com lm)

2. testes de hipótese

Vamos discutir vários testes de hipóteses usando tanto as funções que implementam estes testes em R quanto simulações. Simulações são uma ferramenta poderosa para testar hipóteses e comparar distribuições.

2.1 Distância à LMC

A distância à Grande Nuvem de Magalhães (LMC) é um “degrau” decisivo na calibração de distâncias. Esta distância é usualmente expressa em termos do “módulo de distância”: \(DM = 5 \log d -5\), onde \(d\) é a distância em parsecs.

O arquivo http://www.explorationday.psu.edu/datasets/LMC_distance.html contém distâncias obtidas por dois conjuntos de métodos diferentes, com indicadores de distância baseados em objetos de população I e de população II.

Será que a média e o desvio padrão de cada conjunto de indicadores são consistentes entre si?

Os DM para cada população são dados abaixo:

popI = c(18.7,18.55,18.55,18.575,18.4,18.42,18.45,18.59,18.471,18.54,18.54,18.64,18.58)
popII = c(18.45,18.45,18.44,18.3,18.38,18.5,18.55,18.52,18.4,18.45,18.55,18.69)

# tamanho das populações
length(popI)
## [1] 13
length(popII)
## [1] 12
# estatísticas: média e desvio padrão de cada população
c(mean(popI),sd(popI))
## [1] 18.53892308  0.08562073
c(mean(popII),sd(popII))
## [1] 18.47333333  0.09930058
# distribuições de DM
par(mfrow = c(1, 2))
hist(popI, col=rgb(1,0,0,0.5),xlim=c(18.2,18.8), ylim=c(0,6), main='popI', xlab='DM',breaks=10)
hist(popII, col=rgb(1,0,0,0.5),xlim=c(18.2,18.8), ylim=c(0,6), main='popII', xlab='DM',breaks=10)

# histogramas sobrepostos
par(mfrow = c(1, 1))
hist(popII, col=rgb(1,0,0,0.5),xlim=c(18.2,18.8), ylim=c(0,6), main='distância à LMC', xlab='DM')
hist(popI, col=rgb(0,0,1,0.5), add=T)
box()

As médias e desvio padrão para as populações I e II são: \[ DM_I = 18.539 \pm 0.086\] e \[ DM_{II} = 18.473 \pm 0.099 \] Vamos usar o R para comparar os dois conjuntos de dados.

Vamos comparar as médias de DM com a função t.test(). Vamos considerar dois casos: supondo variâncias diferentes e supondo variâncias iguais, para os dois conjuntos de dados.

H\(_0\): as duas populações têm a mesma média

# teste t supondo a mesma variância para as duas populações 
t.test(popI,popII, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  popI and popII
## t = 1.7729, df = 23, p-value = 0.08949
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01094255  0.14212203
## sample estimates:
## mean of x mean of y 
##  18.53892  18.47333
# teste t supondo variâncias diferentes para as duas populações 
t.test(popI,popII, var.equal=FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  popI and popII
## t = 1.762, df = 21.847, p-value = 0.09205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0116397  0.1428192
## sample estimates:
## mean of x mean of y 
##  18.53892  18.47333

O valor p, nos dois casos, é relativamente alto (\(\sim10\)%). Não se pode rejeitar H\(_0\) com \(\alpha=0.01\) mas se pode rejeitar com \(\alpha=0.1\).

Note que no primeiro caso foi aplicada a versão de Student do teste e, no segundo, a versão de Welch: a versão de Student assume distribuições gaussianas de mesma variância, enquanto que na versão de Welch as variâncias podem ser diferentes.

Vamos agora comparar as variâncias com o teste F, usando a função var.test().

H\(_0\): a variância das duas amostras é igual

var.test(popI,popII)
## 
##  F test to compare two variances
## 
## data:  popI and popII
## F = 0.74345, num df = 12, denom df = 11, p-value = 0.6169
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.216775 2.469370
## sample estimates:
## ratio of variances 
##          0.7434543

o valor p é muito alto: não se pode rejeitar a hipótese de que as duas variâncias são iguais.

2.2 Testes com simulações de bootstrap

O bootstrap é uma técnica poderosa de simulação de dados e pode ser usado para se fazer testes estatísticos.

Por exemplo, suponha que queremos testar a hipótese de que a mediana dos DM da população I é maior que a dos da população II:

# medianas de DM observadas
c(median(popI),median(popII))
## [1] 18.55 18.45

Vamos fazer isso com simulações:

# reprodutibilidade
set.seed(123)

# tamanho das amostras
nI = length(popI)
nII = length(popII)

# número de simulações
nsim = 10000

# simulação 
simI = replicate(nsim, median(sample(popI,size=nI, replace=TRUE)))


simII = replicate(nsim, median(sample(popII,size=nII, replace=TRUE)))


par(mfrow = c(1, 2))
hist(simI,col='aquamarine',xlab='mediana de DM', main='pop I')
hist(simII,col='aquamarine',xlab='mediana de DM', main='pop II')

sim = simI - simII
par(mfrow = c(1, 1))
hist(sim,col='orange',xlab='diferença nas medianas', main='med(pop I) - med(popII)')

# probabilidade de que a mediana de DM da população I é maior que a da população II:
prob = length(sim[sim > 0]) / nsim

prob
## [1] 0.9545

A probabilidade da mediana dos indicadores de população I ser maior que a dos de população II é de 95.5%.

2.3 Einstein x Newton

Se a variável tem uma distribuição mais ou menos gaussiana, pode-se usar o teste-t (na verdade se usa ele mesmo sem essa condição!).

O desvio gravitacional da luz de uma estrela pelo Sol, previsto pela Relatividade Geral, é 1.75 arcsec. Durante um eclipse um astrônomo mede 3 estrelas e obtém \(1.60 \pm 0.10\) arcsec. Esse resultado é consistente com a hipótese nula de que \(\alpha=1.75\) arcsec?

# vamos calcular a estatística t:
alfa_GR = 1.75
alfa_obs = 1.60
sigma = 0.10
n = 3

t = (alfa_obs-alfa_GR)/(sigma/sqrt(n))
t
## [1] -2.598076
# valor-p:
1-pt(t,df=n-1)
## [1] 0.939155

O valor-p é alto e não permite descartar a Relatividade Geral

E a teoria newtoniana? Nesse caso o valor previsto para a deflexão da luz é exatamente a metade do valor da relatividade geral:

alfa_N = 1.75/2
alfa_N
## [1] 0.875
# vamos calcular a estatística t:
t = (alfa_obs-alfa_N)/(sigma/sqrt(n))

# valor-p:
1-pt(t,df=n-1)
## [1] 0.003140981

Já este valor de p é bem mais baixo! As observações não são consistentes com a teoria da gravitação newtoniana com um nível de significância de 1%.

2.4 Teste de proporção numa amostra

Vamos supor que selecionamos ao acaso 100 objetos de uma amostra de estrelas e galáxias e encontramos entre eles 42 galáxias. Esse valor é consistente com a proporção verdadeira \(\theta\) ser 50%?

Hipótese nula H\(_0\): \(\theta=0.5\)

Vamos fazer um teste de proporção via simulação

Nesse caso, simula-se distribuições de acordo com \(H_0\) e verifica-se a probabilidade de se obter uma proporção \(\theta\)

Assume-se que as distribuições são binomiais.

Vamos exemplificar com o caso de 42 galáxias em uma amostra de 100 objetos.

# baseado em
# http://mlwiki.org/index.php/Simulation_For_Proportions

# reprodutibilidade
set.seed(123)

# H0: proporção de galáxias
theta=0.5

# tamanho da amostra
n=100

# número de simulações
nsim = 100000

# número de galáxias observado
nobs = 42

# proporção observada
theta_obs = nobs/n

# simulação
theta_sim = rbinom(n=nsim, size=n, prob=theta) / n

# valor-p
valor.p = sum(theta_sim <= theta_obs) / nsim

valor.p
## [1] 0.06651

Vamos comparar o valor-p obtido com a simulação com o da função prop.test():

prop.test(42,100,p=theta,alternative='less')
## 
##  1-sample proportions test with continuity correction
## 
## data:  42 out of 100, null probability theta
## X-squared = 2.25, df = 1, p-value = 0.06681
## alternative hypothesis: true p is less than 0.5
## 95 percent confidence interval:
##  0.0000000 0.5072341
## sample estimates:
##    p 
## 0.42
# o valor p dado pelo teste é muito parecido com o obtido pela simulação

2.5 Exercícios com a biblioteca MASS do R

Visitem o site http://www.r-tutor.com/category/r-packages/mass

ele contém vários tutoriais mostrando como aplicar testes estatísticos com a biblioteca MASS:

  • Wilcoxon Signed-Rank Test: teste não paramétrico para testar se amostras vêm da mesma população
  • teste de independência de duas variáveis
  • proporção de duas populações
  • etc…

3 comparação de distribuições

Os testes vão depender do tipo das distribuições

3.1 teste do \(\chi^2\) para comparar variáveis nominais

Novais & Sodré (2018) classificaram galáxias do projeto CALIFA de acordo com a morfologia da emissão em H\(\alpha\) em duas grandes classes: C (de central) e Ex (de extended), quando o máximo da emissão ocorre no centro da galáxia e fora dele, respectivamente.

Vamos examinar se esses tipos se correlacionam de algum modo com os tipos morfológicos das galáxias. A tabela abaixo mostra a frequência das classes na amostra:

C Ex
E-S0 29 1
S-early 12 11
S-late 22 11

Duas variáveis aleatórias x e y são ditas independentes se a probabilidade de uma não é afetada pela outra.

Hipótese nula: \(H_0\): as duas variáveis (tipo morfológico e tipo de perfil H\(\alpha\)) são independentes

Hipótese alternativa: \(H_A\): as duas variáveis são relacionadas

Vamos construir uma matriz desses dados como

NIJ = matrix(c(29,12,22,1,11,11),nrow=3)
NIJ
##      [,1] [,2]
## [1,]   29    1
## [2,]   12   11
## [3,]   22   11

estatística: \[ \chi^2 = \sum_{i,j} {(N_{ij} -n_{ij})^2 \over n_{ij}} \] onde \(N_{ij}\) são as contagens observadas e \(n_{ij}\) é o valor esperado para as contagens. Se as linhas e colunas forem independentes pode-se mostrar que: \[n_{ij} = {(\sum_j N_{ij}) (\sum_i N_{ij}) \over N} \]

Assim,

ni = apply(NIJ,1,sum)
nj = apply(NIJ,2,sum)
chi2=0
for (i in 1:3){
  for (j in 1:2){
    nij=ni[i]*nj[j]/sum(NIJ)
    chi2=chi2+(NIJ[i,j]-nij)^2/nij
  }
}

#chi2
chi2
## [1] 14.34133
# número de graus de liberdade neste caso:
ndof=length(ni)*length(nj)-length(ni)-length(nj)+1
ndof
## [1] 2
# chi2 reduzido
chi2/ndof
## [1] 7.170666
# valor p
1-pchisq(chi2,ndof)
## [1] 0.0007688102

ou podemos usar diretamente a função chisq.test() do R:

chisq.test(NIJ)
## 
##  Pearson's Chi-squared test
## 
## data:  NIJ
## X-squared = 14.341, df = 2, p-value = 0.0007688

Como o valor p é pequeno, (menor que um nível de significância de, digamos, \(\alpha = 0.01\)), rejeitamos \(H_0\) e concluímos que o teste indica que as duas variáveis são relacionadas.

3.2 teste de Kolmogorov-Smirnov (KS)

KS é um teste muito usado para verificar se duas distribuições unidimensionais são provenientes de uma mesma distribuição de probabilidades (a hipótese nula), comparando suas distribuições cumulativas.

Vamos ilustrar os indicadores de distância baseados em objetos de população I e de população 2 para a LMC:

ks.test(popI,popII, alternative="two.sided")
## Warning in ks.test(popI, popII, alternative = "two.sided"): cannot compute
## exact p-value with ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  popI and popII
## D = 0.44231, p-value = 0.1739
## alternative hypothesis: two-sided

Como o teste KS é baseado no ordenamento dos valores das duas distribuições, ele avisa quando há dados repetidos (ties) e usa uma aproximação para estimar o valor-p.

No caso, o valor p é alto: não se pode rejeitar \(H_0\), as duas distribuições são consistentes entre si!

A função ks.test() também faz o chamado one sample test, onde se pode verificar se a amostra em análise é consistente com uma certa distribuição, usando sua forma cumulativa. Por exemplo, será que a distribuição dos módulos de distância dos objetos de população I é gaussiana?

hist(popI,col='gold')

ks.test(popI,'pnorm')
## Warning in ks.test(popI, "pnorm"): ties should not be present for the
## Kolmogorov-Smirnov test
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  popI
## D = 1, p-value = 1.022e-11
## alternative hypothesis: two-sided

O valor-p é muito pequeno: os dados não são consistentes com uma gaussiana.

4. análise de correlação

Vamos considerar o uso dos coeficientes de correlação para analisar a relação entre duas variáveis.

Funções úteis:

  • cor(): correlação entre dois vetores, todas as colunas de um data.frame ou dois data.frames
  • cor.test(): testa a significância de uma correlação; a hipótese nula é que não há correlação

Vamos ilustrar estas funções com dados de parâmetros de galáxias em Novais & Sodré (2018). Será que a idade média e a metalicidade das populações estelares das galáxias da amostra estão correlacionadas?

# arquivo de dados
# coloque aqui o endereço onde está a tabela
tab = read.table('/home/laerte/cursos/dados/2022/dados/novais.txt', head = TRUE)
head(tab)
##    Galaxy     Mr  u.r logt  logZ logMstar Hubble_type Morphological_class
## 1  IC0776 -18.69 1.94 8.34 -0.27     9.59          Sd               Slate
## 2  IC1256 -20.81 2.21 9.08 -0.29    10.72          Sb              Searly
## 3  IC1683 -20.75 2.54 9.31 -0.24    10.76          Sb              Searly
## 4  IC4566 -21.51 2.88 9.30 -0.26    11.01          Sb              Searly
## 5 NGC0001 -21.11 2.24 8.89 -0.35    10.82         Sbc               Slate
## 6 NGC0036 -21.86 2.48 9.30 -0.34    11.22          Sb              Searly
##   Halpha_profile   cr ceHalpha ccHalpha  eps
## 1             CL 1.95     0.55     0.20 0.28
## 2             CL 2.19     0.57     0.22 0.20
## 3             CL 2.59     0.30     0.77 0.28
## 4             EX 2.24     0.63     0.08 0.26
## 5             CE 3.04     0.52     0.40 0.25
## 6             EX 2.49     0.72     0.05 0.25
# significado das colunas:
# galaxy name, absolute magnitude $M_r$, colour u−r, mean age (log(age/yr)), metallicity (Z/Zsun), stellar mass (log(M/Msun)), Hubble type, morphological class, type of the Halpha profile, light concentration, effective concentration, central concentration, and ellipticity

# a correlação entre dois vetores pode ser calculada com a função cor()
# faça ?cor para saber mais sobre esta função

# o parâmetro method permite escolher a estatística de correlação e covariância, onde as opções são "pearson" (default), "spearman" ou "kendall"
cor(tab$logt,tab$logZ,method = 'pearson')
## [1] 0.4555887
cor.test(tab$logt,tab$logZ,method = 'pearson')
## 
##  Pearson's product-moment correlation
## 
## data:  tab$logt and tab$logZ
## t = 4.6906, df = 84, p-value = 1.045e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2697521 0.6087068
## sample estimates:
##       cor 
## 0.4555887
cor.test(tab$logt,tab$logZ,method = 'spearman')
## Warning in cor.test.default(tab$logt, tab$logZ, method = "spearman"): Cannot
## compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  tab$logt and tab$logZ
## S = 50877, p-value = 2.879e-07
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5200067
cor.test(tab$logt,tab$logZ,method = 'kendall')
## 
##  Kendall's rank correlation tau
## 
## data:  tab$logt and tab$logZ
## z = 4.7678, p-value = 1.862e-06
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.3568582
plot(tab$logt,tab$logZ,xlab='log t (Ganos)',ylab='log Z/Zsol',main='',pch=20,col='salmon')

Note como diferentes estatísticas produzem resultados diferentes para o coeficiente de correlação.

No caso, existe correlação, ela não é muito forte (os coeficientes de correlação não são muito altos) mas é significativa (valor-p pequeno).

Vamos agora avaliar a correlação entre a magnitude absoluta na banda r e a massa estelar:

cor.test(tab$Mr,tab$logMstar,method = 'spearman')
## Warning in cor.test.default(tab$Mr, tab$logMstar, method = "spearman"): Cannot
## compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  tab$Mr and tab$logMstar
## S = 210120, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.9823558
plot(tab$Mr,tab$logMstar,xlab='Mr',ylab='log Mstar/Msol',main='',pch=20,col='salmon')

essa (anti)correlação é muito forte!

Pode-se calcular várias correlações ao mesmo tempo:

library(corrplot)
## corrplot 0.92 loaded
library(RColorBrewer)
M <-cor(tab[,c(2:6)],method = 'spearman')  
M
##                  Mr        u.r       logt       logZ   logMstar
## Mr        1.0000000 -0.5844048 -0.6181165 -0.4360995 -0.9823558
## u.r      -0.5844048  1.0000000  0.8135370  0.4571195  0.6258605
## logt     -0.6181165  0.8135370  1.0000000  0.5200067  0.6826351
## logZ     -0.4360995  0.4571195  0.5200067  1.0000000  0.4785962
## logMstar -0.9823558  0.6258605  0.6826351  0.4785962  1.0000000

A função corrplot da biblioteca de mesmo nome permite plotar a significância de várias variáveis ao mesmo tempo:

library(corrplot)
library(RColorBrewer)
M <-cor(tab[,c(2:6,10:13)],method = 'pearson')  
corrplot(M, type="upper", order="hclust",
         col=brewer.pal(n=8, name="RdYlBu"))

5. o caso da galáxia sem matéria escura

van Dokkum et al. (2018), Nature: A galaxy lacking dark matter

van Dokkum e colaboradores publicaram em 2018 um artigo na revista Nature onde argumentavam que a galáxia ultra difusa NGC1052-DF2 não apresenta matéria escura.

Para chegar a essa conclusão os autores mediram a velocidade radial de 10 aglomerados globulares dessa galáxia. A dispersão dessas velocidades junto com medidas das posições dos aglomerados permite estimar a massa da galáxia fazendo uso do teorema do virial.

O teorema do virial estabelece que, em um sistema gravitacional em equilíbrio, a energia cinética é metade do módulo da energia potencial: \[T = {1 \over 2} M <v^2> ~=~ {1 \over 2} |W|~=~ {G M^2 \over 2 R}~~~~\longrightarrow ~~~~ M = {R <v^2> \over G }\] Assim, fazendo-se medidas de velocidades e posições, pode-se determinar a massa do sistema. O termo \(<v^2>\) é a dispersão de velocidades. Foi aplicando o teorema do virial nos aglomerados de galáxias de Virgo e Coma que Zwicky, em 1933, descobriu a matéria escura.

No paper os autores escrevem: “We infer that its velocity dispersion is less than 10.5 kilometres per second with 90 per cent confidence” e é isso que acaba levando à conclusão de que a galáxia não tem matéria escura.

O resultado chamou muito a atenção e, em particular, Nicolas Martin escreveu no twiter (https://twitter.com/nfmartin1980/status/982245161735372804)

Como escreveu o Nicolas Martin, “extraordinary claims requires extraordinary proofs”, e ele e um aluno reexaminaram o que os autores fizeram. Em resumo: a análise estatística do artigo está errada. Aparentemente os autores não levaram em conta os grandes erros nas medidas das velocidades (Nature, que papelão!)… O Nicolas enfatiza a necessidade de sempre se deixar claro o modelo estatístico que se usa! E, como veremos, não se deve aplicar cegamente testes estatísticos.

Os dados de velocidade dos aglomerados globulares (em relação ao centro da galáxia) e respectivos erros são

v=c(15,-4,2,11,1,2,-1,-14,-39,-3)
err=c(7,15,7,3,6,5,10,6,12,13)
#o histograma dessas velocidades e de seus erros  é
par(mfrow=c(1,2))
hist(v,col='aquamarine',main='v')
hist(err,col='aquamarine',main='erros em v')

# e a média e a dispersão de v são:
mean(v)
## [1] -3
sd(v)
## [1] 14.9369

Notem que os erros são, em alguns casos, da mesma ordem de grandeza que as velocidades! Medir com precisão essas velocidades é difícil!

no paper: “We infer that its velocity dispersion is less than 10.5 kilometres per second with 90 per cent confidence”

A pergunta que os autores do paper provavelmente fizeram foi algo do tipo: qual é a probabilidade de se ter um valor para a dispersão de velocidades igual ou menor que 10.5 km/s?

Note que, como a massa é proporcional à dispersão de velocidades, um limite em \(<v^2>\) é um limite em massa.

Supondo que o ponto mais a esquerda é um outlier e removendo-o da análise, temos que

v=c(15,-4,2,11,1,2,-1,-14,-3)
err=c(7,15,7,3,6,5,10,6,13)

#o histograma dessas velocidades é
hist(v,col='aquamarine',main='velocidades')

# e a média e a dispersão de v são:
mean(v)
## [1] 1
sd(v)
## [1] 8.42615

Vamos inicialmente supor que os dados foram extraídos de uma gaussiana de média 1 km/s e desvio padrão 8.4 km/s e vamos determinar qual é a probabilidade de se obter uma distribuição com dispersão de velocidades menor que 10.5 km/s fazendo uma simulação:

# reprodutibilidade
set.seed(123)

# número de simulações
nsim = 100000

# parâmetros da distribuição de velocidades
n = length(v)
vmed = mean(v)
sdvm=sd(v)

# velocidade limite
sigvlim=10.5

# simulação
sigv = replicate(nsim,sd(rnorm(n, mean=vmed,sd=sdvm)))

hist(sigv,col='aquamarine',main='sigma_v')
abline(v=sigvlim,lwd=2)

# probabilidade:
length(sigv[sigv < sigvlim])/nsim
## [1] 0.86592

Não foi bem isso que os autores fizeram mas o resultado é parecido com o que eles obtiveram.

Vamos agora levar em consideração o erro nas velocidades e supor que cada velocidade é produzida por uma distribuição N(mean=v,sd=err). Vamos gerar muitas realizações dessas medidas, determinar a dispersão de velocidades em cada caso e analisar a distribuição resultante das dispersões de velocidade.

Primeiro sem o outlier:

# reprodutibilidade
set.seed(123)

# número de simulações
nsim = 100000

# simulação
sigv = replicate(nsim,sd(rnorm(n, mean=v,sd=err)))

hist(sigv,col='aquamarine',main='sigma_v')
abline(v=sigvlim,lwd=2)

# probabilidade:
length(sigv[sigv < sigvlim])/nsim
## [1] 0.30771

Agora incluindo o “outlier”:

vm=c(15,-4,2,11,1,2,-1,-14,-39,-3)
errv=c(7,15,7,3,6,5,10,6,12,13)


ns = length(vm)

e refazendo a simulação:

# reprodutibilidade
set.seed(123)

# número de simulações
nsim = 100000

# simulação
sigv = replicate(nsim,sd(rnorm(ns, mean=vm,sd=errv)))

hist(sigv,col='aquamarine',main='sigma_v')
abline(v=sigvlim,lwd=2)

# probabilidade:
length(sigv[sigv < sigvlim])/nsim
## [1] 0.01297

Com ou sem outlier as probabilidades são muito menores que 90%!

Não se deve aplicar testes estatísticos cegamente!

exercícios

  1. Um estudante precisa testar as hipóteses H\(_0\): \(\mu= 80\) e H\(_A\): \(\mu > 80\) com \(\alpha=0.05\). Analisando a amostra, ele calcula o valor p, 0.214,e conclui que “este resultado prova que H\(_0\) é verdadeiro’’. Comente esta conclusão e a reescreva corretamente.
  1. Observamos algumas estrelas que achamos que podem ser variáveis e uma outra de comparação, que achamos que não é. Observamos 10 noites seguidas obtendo as seguintes medidas (os vetores correspondem a noites sucessivas):
estrela1 = c(14.99, 14.73, 15.18, 15.10, 15.20, 14.89, 15.50, 14.69, 15.25, 14.85)  
estrela2 = c(15.55, 15.48, 15.68, 15.53, 15.43, 15.54, 15.70, 15.57, 15.60, 15.65) 
comparacao = c(15.28, 15.30, 15.32, 15.35, 15.26, 15.28, 15.30, 15.27, 15.26, 15.28)
  1. Faça box plots para visualizar as distribuições das magnitudes de cada um desses objetos.
  1. Adotando um nível de confiança de 99%, o que você concluiria sobre a variabilidade dessas duas estrelas? Use tanto a função var.test() quanto simulações e comente seu resultado
  1. A função ks.test do R permite fazer testes comparando duas amostras (two-sample test) ou comparando os dados com uma distribuição teórica (one-sample test). Considere a sequência de pontos: 1.41, 0.25, 0.26, 1.97, 0.33, 0.55, 0.77, 1.46, 1.18, 0.23. Discuta se existe alguma evidência de que estes dados não resultem de uma distribuição uniforme entre 0 e 2. Faça um teste de Kolmogorov-Smirnov e comente o resultado.
  1. Examine a correlação entre a cor (u-r) com a idade (logt) e com a metalicidade (logZ) para a amostra em novais.txt. O que você conclui? Discuta se e como esse resultado depende do método usado para medir a correlação.