- regressão linear
- relações de escala de galáxias early-type
- otimização de modelos não-lineares
- testes de hipóteses
- comparação de distribuições
- exercícios
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)
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.056106Vamos 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))
\[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)
text(10.5,0.4,'início')
text(3,1.7,'fim')
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
É 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 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)
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))
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áxiasestas 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\).
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 é 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/).
# 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.
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
NULLvamos 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.06764456sd(p2)[1] 9.183396sd(p3)[1] 0.1098924Compare 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 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.
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] 13length(popII)[1] 12# estatísticas: média e desvio padrão de cada população
c(mean(popI),sd(popI))[1] 18.53892308 0.08562073c(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.
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.45Vamos 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.9545A probabilidade da mediana dos indicadores de população I ser maior que a dos de população II é de 95.5%.
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…
Os testes vão depender do tipo das distribuições
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 11estatí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.0007688102ou 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.0007688Como 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.
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-sidedComo 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-sidedO valor-p é muito pequeno: os dados não são consistentes com uma gaussiana.
- 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.
- 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)
- Faça box plots para visualizar as distribuições das magnitudes de cada um desses objetos.
- 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().