aula de hoje:

    1. regressão linear
    1. relações de escala de galáxias early-type
    1. otimização de modelos não-lineares
    1. testes de hipóteses
    1. comparação de distribuições
    1. análise de correlação
    1. exercícios

1. regressão linear

regressão linear \(=\) ajuste de uma reta

Vamos gerar uma amostra de 10 pontos \((x,y)\) obedecendo uma relação \(y=a+bx\), com \(a=1\) e \(b=2\) mais ruído gaussiano com desvio padrão 1.8 e com \(0 < x < 15\):

# reprodutibilidade
set.seed(123)

# 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

# incluindo erros
err=rnorm(n,mean=0,sd=1.8)
y=y+err

# visualização
plot(x,y,main='y=1+2x     sd=1.8',pch=20, col='red',cex=2)

# modelo da reta:
abline(a=1,b=2,lwd=2)

estimativa dos parâmetros de um modelo linear:

Dados \(n\) pontos \((x_i,y_i)\) com erros (gaussianos) \(\sigma_i\), a verossimilhança é \[log {\cal L}(a,b) \propto -{1 \over 2}\chi^2,\] onde \[\chi^2 (a,b) = \sum_{i=1}^n {(y_i-a-b x_i)^2 \over \sigma_i^2}.\]

A solução de máxima verossimilhança pode ser determinada analiticamente: \[ a = {S_{xx} S_y - S_x S_{xy} \over \Delta } ~~~~~~~~~~~ b = {S S_{xy}-S_x S_y \over \Delta} \] onde \[ S = \sum_{i=1}^N {1 \over \sigma_i^2} ~~~~~~~~~ S_x = \sum_{i=1}^N {x_i \over \sigma_i^2} ~~~~~~~~~ S_y = \sum_{i=1}^N {y_i \over \sigma_i^2} ~~~~~~~~~ S_{xx} = \sum_{i=1}^N {x_i^2 \over \sigma_i^2} ~~~~~~~~~ S_{xy} = \sum_{i=1}^N {x_i y_i \over \sigma_i^2}\] \[ \Delta = S S_{xx} - S_x^2 \]

Os erros nos parâmetros \(a\) e \(b\) são: \[ \sigma_a^2 ={S_{xx} \over \Delta}~~~~~~~~~~~~~~~~~ \sigma_b^2 ={S \over \Delta} \]

S=sum(1/err^2)
Sx=sum(x/err^2)
Sy=sum(y/err^2)
Sxx=sum(x^2/err^2)
Sxy=sum(x*y/err^2)

Delta=S*Sxx-Sx^2
a=(Sxx*Sy-Sx*Sxy)/Delta
b=(S*Sxy-Sx*Sy)/Delta
siga=sqrt(Sxx/Delta)
sigb=sqrt(S/Delta)

#parâmetros ajustados
#comparem-os com os usados na simulação dos dados: a diferença é devida ao ruído.
a; b
## [1] 1.635979
## [1] 1.949556
#erros nos parâmetros
siga; sigb
## [1] 0.7743093
## [1] 0.0846401
# note que o erro na inclinação é muito menor que no "ponto zero"

#visualização
xp=seq(0,15,0.1)
yp=a+b*xp
plot(x,y,main='modelo: y=1+2x',col='red',pch=20,cex=2,lwd=2,ylim=c(0,30))

#linha verde: reta ajustada
lines(xp,yp,col='green')

#linha preta: modelo inicial
yp=afit+bfit*xp
lines(xp,yp,col='black',lwd=2)

legend('topleft',c('ajuste','modelo inicial'),col=c('green','black'),pch = c(20,20))

# a reta simulada e a inicial ficam praticamente uma em cima da outra nesse caso

# chi2
chi2=sum((y-a-b*x)^2/err^2)
chi2
## [1] 8.448851
# chi2 reduzido
# nu: número de graus de liberdade = número de pontos menos o número de parâmetros do modelo
nu=n-2
chi2/nu
## [1] 1.056106

Podemos também resolver este problema usando matrizes

(método de Gauss-Markov):

Vamos escrever as \(n\) equações do ajuste linear como matrizes: \[ \begin{bmatrix} y_1 \\ y_2 \\ ... \\ y_n \end{bmatrix} = \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ ... \\ 1 & x_n \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} \] ou, de forma mais compacta, \[\vec{y} = \mathbf{M} \vec{w},\] onde \(\vec{y}\) é o vetor da medidas \(y\), \(\vec{w}\) é o vetor de parâmetros.

Vamos multiplicar ambos os lados de \(\vec{y} = \mathbf{M} \vec{w}\) pela matriz transposta de \(\mathbf{M}\): \[\mathbf{M}^T \vec{y} = \mathbf{M}^T \mathbf{M} \vec{w}\] e em seguida pela matriz inversa de \(\mathbf{M}^T \mathbf{M}\), obtendo \[\vec{w} = (\mathbf{M}^T \mathbf{M})^{-1}\mathbf{M}^T \vec{y}.\] Pode-se mostrar que esta solução é idêntica à de máxima verossimilhança (com erros homocedásticos)!

Note que \[\mathbf{M}^T \vec{y} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \end{bmatrix} \]

\[\mathbf{M}^T \mathbf{M} = \begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} \]

\[(\mathbf{M}^T \mathbf{M})^{-1} = \begin{bmatrix} {\sum x_i^2 \over n\sum x_i^2 -(\sum x_i)^2 } & {-\sum x_i \over n\sum x_i^2 -(\sum x_i)^2 } \\ {-\sum x_i \over n\sum x_i^2 -(\sum x_i)^2 } & {n \over n\sum x_i^2 -(\sum x_i)^2 } \end{bmatrix}\]

A matriz de covariância é \[\mathbf{C} = \sigma^2 (\mathbf{M}^T \mathbf{M})^{-1}\]

Este resultado assume erros homocedásticos de variância \(\sigma^2\). No caso de erros heterocedásticos, pode-se mostrar que \[\vec{w} = (\mathbf{M}^T \Omega \mathbf{M})^{-1}\mathbf{M}^T \Omega \vec{y},\] onde \[\Omega = \begin{bmatrix} 1/\sigma_1^2 & 0 & ... & 0 \\ 0 & 1/\sigma_2^2 & ... & 0 \\ ... & ... & ... & 0 \\ 0 & 0 & ... & 1/\sigma_n^2 \end{bmatrix} \]

Vamos resolver o ajuste linear acima com matrizes:

# solução homocedástica
# matriz M
# cbind: column bind
M = cbind(rep(1,n),x)
# transposta de M
Mt = t(M)
# produto Mt x M
MtM = Mt %*% M
# inversa de MtM
Mi = solve(MtM)
# parâmetros
w = Mi %*% Mt %*% y
# solução homocedástica
w
##       [,1]
##   2.618906
## x 1.840705
# solução heteroscedástica
Omega = diag(n)*1/err^2
wl = solve(Mt %*% (Omega %*% M)) %*% t(M) %*% (Omega %*% y)
wl
##       [,1]
##   1.635979
## x 1.949556
# visualização
plot(x,y,main='modelo: y=1+2x',col='red',pch=20,cex=2,lwd=2,ylim=c(0,30))
yp=w[1]+w[2]*xp
lines(xp,yp,col='blue',lwd=2)
ypl=wl[1]+wl[2]*xp
lines(xp,ypl,col='green',lwd=2)
legend('topleft',c('ajuste heterocedástico','ajuste homocedástico'),col=c('green','blue'),pch = c(20,20))

Vamos agora ilustrar o algoritmo de descida do gradiente com a regressão linear:

\[w_{n+1} = w_{n} - \lambda \nabla \chi^2(w)~~~~~~~~~(0 < \lambda < 1)\]

# parâmetro de aprendizado
lambda <- 0.01

# número de iterações do algoritmo
num_iters <- 2000

# y: variável dependente
# ycalc = X %*% w

# função custo: erro quadrático
custo <- function(X, y, w) sum( (X %*% w - y)^2 ) 

# vetores auxiliares para guardar a história da função de custo e dos parâmetros:
custo_historia <- rep(0,num_iters)
w_historia <- list(num_iters)
ait = NULL
bit = NULL

# inicializa os parâmetros 
w <- matrix(c(10,-10), nrow=2)

# adiciona a coluna de '1's para os valores de x que multiplicam w1
X <- cbind(1, matrix(x))

# gradiente descendente:
for (i in 1:num_iters) {
  erro <- (X %*% w - y)
  delta <- t(X) %*% erro / length(y)
  w <- w - lambda * delta
  custo_historia[i] <- custo(X, y, w)
  w_historia[[i]] <- w
  ait[i]=w[1]
  bit[i]=w[2]
}

# resultado
print(w)
##          [,1]
## [1,] 2.814681
## [2,] 1.822374
# compare os valores dos parâmetros ajustados com os usados para simular os dados

# visualização do resultado:
par(mfrow=c(1,1))
plot(x,y, col=rgb(0.2,0.4,0.6,0.8), main='evolução do ajuste com "descendo o gradiente"')
for (i in c(1,3,6,10,14,seq(20,num_iters,by=10))) {
  abline(coef=w_historia[[i]], col=rgb(0.8,0,0,0.3))
}
abline(coef=w, col='blue')

# trajetória no espaço de parâmetros:
plot(ait,bit,type='l',xlab='a',ylab='b',col='magenta',lwd=3)
points(w[1],w[2],pch=20,cex=2)

A função \(lm()\)

Para ajustar uma reta aos pontos- a regressão linear- usa-se em R a função \(lm\) (de linear model), que tem a sintaxe \[ lm (y \sim model) \] onde \(y\) é um objeto em R contendo a variável dependente e \(model\) a fórmula em R para o modelo matemático adotado.

No caso de uma regressão linear simples, como no nosso exemplo, onde \(x\) é a variável independente e \(y\) a variável dependente, \(y = \beta_0+\beta_1 x\), usa-se \(lm (y \sim x)\).

Importante: lembre-se que ‘linear’ não significa apenas a reta, mas linear nos parâmetros.

Exemplo: ajuste de uma reta aos dados poderados pelo inverso do quadrado dos erros:

lm(y~x,weights = 1/err^2)
## 
## Call:
## lm(formula = y ~ x, weights = 1/err^2)
## 
## Coefficients:
## (Intercept)            x  
##       1.636        1.950

Faça um ?lm para saber mais sobre essa função e sua sintaxe

resultado da regressão:

É melhor se atribuir um objeto a \(lm()\) pois ele pode ser usado com outros comandos:

# crio o objeto ajuste
ajuste=lm(y~x,weights = 1/err^2)
ajuste
## 
## Call:
## lm(formula = y ~ x, weights = 1/err^2)
## 
## Coefficients:
## (Intercept)            x  
##       1.636        1.950
#summary
summary(ajuste)
## 
## Call:
## lm(formula = y ~ x, weights = 1/err^2)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2903 -0.9569  0.2686  0.8301  1.0544 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.63598    0.79573   2.056   0.0738 .  
## x            1.94956    0.08698  22.413 1.66e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.028 on 8 degrees of freedom
## Multiple R-squared:  0.9843, Adjusted R-squared:  0.9824 
## F-statistic: 502.4 on 1 and 8 DF,  p-value: 1.661e-08
#significado do resultado do summary do lm:
# 1. fórmula usada no cálculo

# 2. summary dos resíduos (observado - calculado) ponderados, se for o caso
###  em um bom ajuste deve haver uma simetria em relação à média ou mediana

# 3. para cada coeficiente: 
###  estimativa
###  erro padrão - média e desvio padrão da distribuição do coeficiente
###  valor t: coeficiente dividido pelo desvio padrão - quanto maior (em módulo) mais significativo o parâmetro é
###  Pr(>|t|): valor p [2*pt(t, df = df.residual(ajuste), lower.tail = FALSE)] - probabilidade de se observar um valor maior ou igual a |t| - quanto menor p, mais improvável que o resultado ocorra por acaso
###  asteriscos: 3 significa que o valor é significativo

# 4. erro padrão dos resíduos: medida da qualidade do ajuste
###  desvio padrão da distribuição dos resíduos

# 5. Multiple R-squared, Adjusted R-squared: medida da qualidade do ajuste- fração da variância explicada

# 6. F-statistics: razão entre a variância explicada pelo modelo (soma dos quadrados da regressão) e a não explicada (soma dos quadrados dos erros)
###  indicador de uma relação entre x e y: quanto maior melhor
#atributos do ajuste:
attributes(ajuste)
## $names
##  [1] "coefficients"  "residuals"     "fitted.values" "effects"      
##  [5] "weights"       "rank"          "assign"        "qr"           
##  [9] "df.residual"   "xlevels"       "call"          "terms"        
## [13] "model"        
## 
## $class
## [1] "lm"
# os atributos permitem obter várias quantidades interessantes do ajuste:
# coeficientes ajustados
a1=ajuste$coefficients[1]
b1=ajuste$coefficients[2]
a1; b1
## (Intercept) 
##    1.635979
##        x 
## 1.949556
# essas quantidades podem ser obtidas também com algumas funções:
# coeficientes ajustados:
coefficients(ajuste)
## (Intercept)           x 
##    1.635979    1.949556
# resíduos do ajuste:
residuals(ajuste)
##           1           2           3           4           5           6 
##  2.66873421  0.79014370 -2.60363571 -1.20417624 -0.72656300  1.60183842 
##           7           8           9          10 
##  0.41127845  0.76066129 -0.01950519 -1.29099421
# valores ajustados para cada x:
fitted.values(ajuste)
##         1         2         3         4         5         6         7         8 
## 10.045708 24.688660 13.595833 27.458363 29.138390  2.968204 17.079551 27.733299 
##         9        10 
## 17.761785 14.988922

calculando \(y\) para novos valores de \(x\) com predict()

Muitas vezes queremos usar o resultado do ajuste para prever valores de \(y\) para novos valores de \(x\).

Isso é feito com a função predict():

# vamos estimar y para os 3 valores de x do vetor xnovo:
xnovo=c(3,5,7.5)

# usa-se a funcão predict():
ynovo=predict(ajuste,data.frame(x = xnovo))

# visualização
plot(x,y,main='y=a+bx=1+2x',pch=20,cex=2,col='red')
abline(ajuste,col='blue',lwd=2)
points(xnovo,ynovo, col="black",pch=20,cex=2)

outros usos de lm()

A figura abaixo ilustra a sintaxe (pouco intuitiva) do parâmetro \(model\) para vários modelos matemáticos. Note que os parâmetros são denominados \(\beta\) e as variáveis independentes A, B, C:

lm() permite implementar vários modelos interessantes!

Vamos ilustrar com o ajuste de uma polinômio de segundo grau (parábola) mais ruído gaussiano:

# reprodutibilidade
set.seed(123)
# número de pontos
N=20
# número de parâmetros
M=3
# parâmetros da simulação
a1=1
a2=-1
a3=0.5
sigma = 0.5
# simulação dos valores de x entre 0 e 5
x=runif(N,min=0,max=5)
# vamos supor que os pontos têm erros gaussianos, gerados com o mesmo desvio padrão sigma
sig=rep(sigma,N)
# simulação da variável dependente
y=a1+a2*x+a3*x^2+rnorm(N,0,sig)
# ajuste dos parâmetros com lm()
parabola = lm(y ~ x +I(x^2))
coefficients(parabola)
## (Intercept)           x      I(x^2) 
##   1.3721938  -1.2047424   0.5138822
# visualização 
plot(x,y,pch=20,cex=2)
xm=seq(0,5,0.1)
ym=a1+a2*xm+a3*xm^2
lines(xm,ym,col='red',lwd=3)
yp=parabola$coefficients[1]+parabola$coefficients[2]*xm+parabola$coefficients[3]*xm^2
lines(xm,yp,col='blue',lwd=3)
legend('topleft',c('modelo','ajuste'),col=c('red','blue'),pch = c(20,20))

2. Relações de escala de galáxias early-type

Galáxias E e S0 obedecem a relações entre seus parâmetros físicos que denominamos relações de escala e que são úteis na determinação da distância a estes sistemas de forma independente do redshift!

Vamos examinar algumas dessas relações usando os dados de galáxias do aglomerado de Virgo que constam da Tabela 2 de Dressler et al. (1987, ApJ, 313,42).

Note que as galáxias em um aglomerado estão todas a mais ou menos a mesma distância de nós.

# leitura de dados
dados <- read.table(file="tab_virgo.dat", header=TRUE)

# dados:
head(dados)
##   ID  NAME   S log_Dn    BT log_Ae Sigmae log_sigma   Mg2
## 1 VI N4239 245  1.258 13.56  1.528  21.69     1.716 0.148
## 2 V2 N4365 254  1.868 10.63  2.068  21.46     2.412 0.311
## 3 V3 N4374 255  2.028 10.15  2.028  20.78     2.480 0.310
## 4 V4 N4387 257  1.478 12.86  1.498  20.84     2.059 0.236
## 5 V5 N4406 258  1.978  9.88  2.278  21.76     2.355 0.299
## 6 V6 N4434 260  1.448 12.80  1.588  21.23     2.009 0.241
# ID - identificacao da galaxia
# NAME - nome da galaxia
# S - outra identificação
# log_Dn - log do diâmetro da abertura circular (em arcsec) contendo o brilho superficial integrado mu_b = 20.75 mag/arcsec^2
# BT - magnitude total (extrapolada) na banda B
# log_ Ae - log do diâmetro (em arcsec) da abertura circular contendo metade da luz da galáxia na banda B - Ae/2 é o raio efetivo
# Sigmae - brilho superficial em mag/arcsec^2 dentro de Ae
# log_sigma - log da dispersão central de velocidades (km/s) da galáxia
# Mg2 - índice espectral Mg2 - indicador de metalicidade, medido no espectro da galáxias

estas galáxias são todas do Aglomerado de Virgo: estão essencialmente à mesma distância da Terra, então as magnitudes absolutas são iguais às aparentes a menos de uma constante (o módulo de distância), que é a mesma para todas as galáxias. O módulo de distância (ou distance modulus) é \[\mu = m-M\] onde \(m\) é a magnitude aparente (uma medida de fluxo) em uma dada banda e \(M\) é a magnitude absoluta (uma medida de luminosidade) na mesma banda

relação de Faber-Jackson: relação entre a luminosidade (ou magnitude absoluta) e a dispersão de velocidades: \[M = a + b \times log\sigma\]

NB: a dispersão de velocidades e a luminosidade de uma galáxia não dependem da distância, mas a magnitude aparente sim; é isso que permite se usar esse tipo de relação (e as demais abaixo) como indicadores de distância

No caso da relação de Faber-Jackson, mede-se a dispersão de velocidades e com ela se determina a luminosidade correspondente; daí, conhecendo-se a magnitude aparente, determina-se a distância.

relação de Kormendy: relação entre o brilho superficial dentro do raio efetivo, \(\Sigma_e\), e o raio efetivo (o raio da abertura circular que contém metade da luminosidade da galáxia) \[log \Sigma_e = a + b \times log R_e\]

como as galáxias estão praticamente à mesma distância de nós, seus diâmetros também serão iguais ao tamanho aparente multiplicados por um fator que é o mesmo para todas as galáxias

a relação de Kormendy se baseia no fato de que, em primeira aproximação, o brilho superficial não depende da distância! Esta relação é muito mais ruidosa que a de Faber-Jackson.

uma relação interessante é o “Plano Fundamental”: um melhor ajuste é obtido combinando-se 3 variáveis: \(BT, \sigma, \Sigma_e\)

Vamos considerar \(BT = f(\Sigma_e,\sigma)\) e usar lm() para fazer um ajuste de mínimos quadrados generalizado:

x = dados$log_sigma
y = dados$Sigmae
z = dados$BT

# mínimos quadrados generalizado
# usa-se uma 'formula' do lm
pf = lm(z ~ x + y)
summary(pf)
## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47029 -0.26853  0.02048  0.20171  0.55111 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.55898    2.13506  17.592 2.41e-12 ***
## x           -5.61719    0.31471 -17.849 1.91e-12 ***
## y           -0.64918    0.09047  -7.175 1.55e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3219 on 17 degrees of freedom
## Multiple R-squared:  0.9518, Adjusted R-squared:  0.9462 
## F-statistic: 167.9 on 2 and 17 DF,  p-value: 6.37e-12
# vizualização dos resultados
par(mfrow=c(1,1))
library("scatterplot3d")
G  <- scatterplot3d(x, y, z, highlight.3d = FALSE, type = "p", main = 'Plano Fundamental',xlab=expression(paste('log ',sigma,' (km/s)')),ylab=expression(paste(Sigma,' (mag/arcsec^2)')),zlab='BT',pch=20)
G$plane3d(pf, draw_polygon = TRUE, draw_lines = FALSE)

suppressPackageStartupMessages(library(FactoClass))

x = dados$log_sigma
y = dados$Sigmae
z = dados$BT
pf = lm(z ~ x + y)
G   <- scatterplot3d(x, y, z, angle = 45, type = "p", main = 'Plano Fundamental',xlab=expression(paste('log ',sigma,' (km/s)')),ylab=expression(paste(Sigma,' (mag/arcsec^2)')),zlab='BT',pch=20)
G$plane3d(pf, lty.box = "solid")
orig <- G$xyz.convert(x, y, z)
plane <- G$xyz.convert(x, y, fitted(pf))
i.negpos <- 1 + (resid(pf) > 0)
segments(orig$x, orig$y, plane$x, plane$y,
         col = c("blue", "red")[i.negpos], 
         lty = 1) # (2:1)[i.negpos]
G <- FactoClass::addgrids3d(x, y, z, angle = 45, grid = c("xy", "xz", "yz"))

as galáxias early-type formam um plano no espaço definido por \(BT\), \(\log \sigma\) e \(\Sigma_e\).

3. Otimização de modelos não-lineares

Vamos agora voltar ao problema específico da minimização do \(\chi^2(\mathbf{w})\) para obtenção dos parâmetros \(\mathbf{w}\) e de seus erros: \[ \chi^2 (\mathbf{w}) = \sum_{i=1}^N {(y_i - f(x_i;\mathbf{w}))^2 \over \sigma_i^2}, \] Vou abordar este problema de duas formas: com a rotina nls e como um problema de otimização, com erros obtidos por bootstrap.

nls: nonlinear least square

  • nls é uma rotina do R para estimativa de parâmetros de modelos não lineares. Sua sintaxe básica é \[ nls(formula, data, start, weights)\]
  • ‘formula’ descreve o modelo não linear, incluindo variáveis e parâmetros, ‘data’ os dados, ‘start’ o ponto inicial e ‘weights’ os pesos (\(1/\sigma_i^2\)), se disponíveis.

Vamos considerar o ajuste do perfil de brilho de uma galáxia elíptica. O brilho superficial \(\mu\) (medido em mag/arcsec\(^2\)) se relaciona com a intensidade \(I\) como \[ \mu (r) = \mu_0 - 2.5 \log I(r)\]

Vamos modelar a intensidade com um perfil de Sérsic: \[ \log I (r) = log I_e + b_n [(r/r_e )^{1/n} - 1] \] onde \(r_e\) é o raio efetivo (o raio que contém metade da luz da galáxia), \(I_e\) é a intensidade em \(r_e\) e \(b_n \simeq 0.868n - 0.142\) para \(0.5 < n < 16.5\).

Discos possuem \(n \sim 1\) (perfis exponenciais), enquanto galáxias elípticas \(n \sim 4\).

O modelo pode ser escrito então como \[ \mu (r) = \mu_0^\prime + 2.5 b_n[(r/r_e )^{1/n} - 1],\] onde \[\mu_0^\prime = \mu_0 -2.5 log I_e.\]

O modelo tem 3 parâmetros: \(\mathbf{w} = \{\mu_0^\prime , n, r_e\}\).

O perfil observado de M49, a galáxia mais brilhante do aglomerado de Virgo, está disponível no site de astroestatística da PennState (https://astrostatistics.psu.edu/).

M49
M49
# leitura dos dados:
NGC4472 <- 
    read.table("NGC4472_profile.dat", header=T)

# vamos olhar a cara dos dados: 
# radius: raio de uma abertura circular em segundos de arco
# surf_mag: brilho superficial em magnitudes por segundos de arco ao quadrado
head(NGC4472)
##   radius surf_mag
## 1  3.526   16.917
## 2  4.015   17.032
## 3  4.557   17.150
## 4  4.966   17.239
## 5  5.491   17.339
## 6  6.039   17.439
# para os objetos na tabela (radius e surf_mag) ficarem acessíveis ao R com os mesmos nomes:
attach(NGC4472)

# ajuste com nls()
NGC4472.fit <-  nls(surf_mag ~ mu0l +2.5*(0.868*n-0.142)*((radius/r.e)^{1/n}-1), data=NGC4472, start=list(mu0l=20.,r.e=120.,n=4.),model=T)

# resultado
summary(NGC4472.fit)
## 
## Formula: surf_mag ~ mu0l + 2.5 * (0.868 * n - 0.142) * ((radius/r.e)^{
##     1/n
## } - 1)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## mu0l  23.35704    0.05478  426.34   <2e-16 ***
## r.e  267.64261    7.05769   37.92   <2e-16 ***
## n      5.95210    0.09718   61.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0377 on 55 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 6.063e-08
# visualização
plot(NGC4472.fit$model$radius, NGC4472.fit$model$surf_mag, pch=20, 
   xlab="r  (arcsec)", ylab=expression(mu ~~ (mag/sq.arcsec)), ylim=c(16,28), 
   cex.lab=1., cex.axis=1., col='red',cex=2)
lines(NGC4472.fit$model$radius, fitted(NGC4472.fit),lwd=2)

atenção: a escolha do ponto inicial (start) pode ser crítica para a convergência do algoritmo!

note que nlm() dá os parâmetros e seus erros.

optim

Vamos agora considerar a minimização do \(\chi^2\) como um problema de otimização.

optim: General-purpose optimization based on Nelder–Mead, quasi-Newton and conjugate-gradient algorithms. It includes an option for box-constrained optimization and simulated annealing

# dimensão dos dados
npt=length(radius)

# função chi2 
# não temos os erros; se os tivéssemos poderíamos incluí-los no cálculo do chi2
fchi2 = function(x) {
s=0
for (k in 1:npt) {
r=radius[k]
bso=surf_mag[k]
bsc=x[1]+2.5*(0.868*x[3]-0.142)*((r/x[2])^{1/x[3]}-1)
s=s+(bso-bsc)^2
}
return(s/2)
}

# chute inicial dos parâmetros
# atenção: a escolha do ponto inicial (start) pode ser crítica para a convergência do algoritmo!
x=c(20,120,4)

# otimização:
res=optim(x,fchi2)
res
## $par
## [1]  23.358306 267.808098   5.954455
## 
## $value
## [1] 0.03908537
## 
## $counts
## function gradient 
##      169       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
# visualização
#plot(raio, bs, pch=20, 
plot(radius, surf_mag, pch=20, 
   xlab="r  (arcsec)", ylab=expression(mu ~~ (mag/sq.arcsec)), ylim=c(16,28), 
   cex.lab=1., cex.axis=1., col='red',cex=2)
   yc = res$par[1]+2.5*(0.868*res$par[3]-0.142)*((radius/res$par[2])^{1/res$par[3]}-1)
lines(radius,yc,col='blue',lwd=2)

Compare os parâmetros obtidos com optim() com os obtidos por nlm()

A função optim determina os parâmetros por otimização de uma função, mas não seus erros. Estes podem ser obtidos por bootstrap. Vamos escrever a função \(\chi^2\) de um jeito útil para depois fazer bootstrap pois, em cada simulação, usa-se uma amostragem diferente dos dados. Para isso vamos definir um vetor \(d\) contendo os índices das observações que entram no cálculo:

# função chi2 escrita de um jeito apropriado para o bootstrap
# d é um vetor com o índice das observações
fchi2 = function(x,d) {
s=0
for (k in 1:npt) {
r=radius[d[k]]
bso=surf_mag[d[k]]
bsc=x[1]+2.5*(0.868*x[3]-0.142)*((r/x[2])^{1/x[3]}-1)
s=s+(bso-bsc)^2
}
return(s/2)
}

# chute inicial dos parâmetros
x=c(10,100,5)

# otimização   d0 é apenas o vetor 1,2...npt indicando que todos os dados da tabela serão usados no chi2
d0=seq(1,npt,1)

# como preciso entrar na optim com um f(x), defino uma nova função
fxyz=function(x){return(fchi2(x,d0))}
# e voilá:
res=optim(x,fxyz)
res
## $par
## [1]  23.357434 267.710532   5.949999
## 
## $value
## [1] 0.03908987
## 
## $counts
## function gradient 
##      254       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

erros por bootstrap

vamos estimar os erros nos parâmetros do ajuste acima com 1000 simulações de bootstrap:

# número de simulações
nboot=1000

# crio uma matriz para armazenar os 3 parâmetros de cada uma das nboot simulações
xb=matrix(rep(0,3*nboot),nboot)

# crio um vetor para armazenar o chi2 minimo de cada bootstrap
chi2min=rep(0,nboot)

# simulações
for (n in 1:nboot) {

# resampling
# 'sample' cria um novo vetor d amostrando d0 com substituição
dboot=sample(d0,npt, replace = TRUE)

# parâmetros iniciais
x=c(20,120,4)

# otimização com os dados sampleados
fxyz=function(x){return(fchi2(x,dboot))}
res=optim(x,fxyz)

# armazena os parâmetros ótimos na matriz xb
for (j in 1:3) {xb[n,j]=res$par[j]}

# armazena o chi2 minimo de cada bootstrap
chi2min[n] = res$value
}
# fim das simulações

# parâmetros das simulações
p1=xb[,1]
p2=xb[,2]
p3=xb[,3]

# desvio padrão dos parâmetros:
sd(p1)
## [1] 0.06764456
sd(p2)
## [1] 9.183396
sd(p3)
## [1] 0.1098924

Compare os erros obtidos por bootstrap com os obtidos por nlm().

Vamos examinar as distribuições dos parâmetros simulados:

par(mfrow=c(2,2))
hist(p1,main='distribuição de Ie',col='red',breaks = 50)
hist(p2,main='distribuição de re',col='blue',breaks = 50)
hist(p3,main='distribuição de n',col='green',breaks = 50)
hist(chi2min,main='distribuição de chi2min', col='yellow',breaks = 50)

# distribuições bi-variadas
par(mfrow=c(2,2))
plot(p1,p2,main='Ie x re',col='red',xlab='mu0',ylab='re',pch=20)
plot(p1,p3,main='Ie x n',col='blue',xlab='mu0',ylab='n',pch=20)
plot(p2,p3,main='re x n',col='green',xlab='re',ylab='n',pch=20)
plot(p3,chi2min,main='n x chi2min',col='salmon',xlab='n',ylab=' chi2min',pch=20)

# como se pode ver, os erros nos parâmetros são correlacionados
# por exemplo:   
cor.test(p1,p2,method='spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  p1 and p2
## S = 1435594, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9913864

4. testes de hipóteses

Vamos discutir alguns 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.

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: para isso se usa o teste F, com 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.

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%.

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…

5. comparação de distribuições

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

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:

tipo 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.

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")
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  popI and popII
## D = 0.44231, p-value = 0.1121
## 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.default(popI, "pnorm"): ties should not be present for the
## one-sample Kolmogorov-Smirnov test
## 
##  Asymptotic 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.

6. 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('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.95 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"))

7. exercícios

  1. A tabela tab_virgo.dat contém também um indicador de metalicidade denominado Mg2, associado a uma banda de absorção do magnésio. Vamos examinar se as propriedades químicas das galáxias early-type, descritas pelo Mg2, podem ser estimadas do brilho superficial ou da dispersão de velocidades. a) Examine a correlação entre Mg2 com o brilho superficial e a dispersão de velocidades. Discuta se e como esse resultado depende do método usado para medir a correlação. b) Use mínimos quadrados para analisar as relações do Mg2 com o brilho superficial e a dispersão de velocidades. Escreva suas conclusões.
  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. Em fotometria diferencial, verifica-se como varia a diferença das magnitudes da candidata a variável e da estrela de comparação. Adotando um nível de confiança de 99%, o que você concluiria sobre a variabilidade das duas estrelas? Use a função var.test().