aula de hoje:

    1. combinação de medidas
    1. propagação dos erros
    1. otimização de funções unidimensionais
    1. otimização de funções bidimensionais
    1. regressão linear
    1. mínimos quadrados generalizado
    1. relações de escala de galáxias early-type
    1. otimização de modelos não-lineares: o perfil de brilho de NGC 4472

1. combinação de medidas

se temos duas medidas de uma certa quantidade e seus erros, \[ x = x_A \pm \sigma_A ~~~~~~~~~~ x = x_B \pm \sigma_B, \] a melhor combinação, pelo método da máxima verossimilhança, é \[ \hat x = {{x_A \over \sigma_A^2} + {x_B \over \sigma_B^2} \over {1 \over \sigma_A^2} + {1 \over \sigma_B^2} } = {w_A x_A + w_B x_B \over w_A + w_B}, \] onde \(w_A = 1/\sigma_A^2\) e \(w_B = 1/\sigma_B^2\)

\(\hat x\) é a média ponderada das medidas, com pesos iguais ao inverso do quadrado dos erros; esta expressão pode ser facilmente generalizada no caso de várias medidas independentes.

para exemplificar, vamos considerar novamente as estimativas de distância à LMC, que constam do arquivo http://www.explorationday.psu.edu/datasets/LMC_distance.html

Temos dois tipos de indicadores de distância, para populações I e II, e queremos combiná-los para obter um módulo de distância ótimo para a LMC.

Os módulos de distância (DM) e seus erros, 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)
errpopI = c(0.16,0.06,0.10,0.20,0.10,0.07,0.07,0.09,0.12,0.10,0.18,0.14,0.05)
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)
errpopII = c(0.09,0.13,0.13,0.14,0.16,0.16,0.19,0.18,0.15,0.16,0.09,0.26)

Notem que, como todas essas medidas são independentes, a melhor estimativa do DM é a média ponderada pelo inverso do quadrados dos erros de todas as medidas individuais, independentemente de serem indicadores de população I ou II:

DM = c(popI,popII)
errDM = c(errpopI,errpopII)

# número de indicadores:
length(DM)
## [1] 25
# média ponderada
DM_MV = weighted.mean(DM,1/errDM^2)
DM_MV 
## [1] 18.50948

Vamos agora ilustrar como utilizar bootstrap para estimar uma barra de erro para esta estimativa.

Bootstrap é uma técnica de estimativa de estatísticas, como variâncias, por reamostragem.

No caso, vamos reamostrar o vetor DM muitas vezes, com substituição, e calcular algumas estatísticas:

# reprodutibilidade:
set.seed(123)

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

# crio um vetor para armazenar o valor médio DM_MV de cada simulação
mediaDM=rep(0,nboot)

# crio um vetor para os indicadores selecionados para cada simulação
ind=seq(1,length(DM),1)

# simulações
for (i in 1:nboot) {
# resampling
# vamos usar sample() para criar um novo vetor de indicadores amostrando ind com substituição
novoind=sample(ind,length(ind), replace = TRUE)

# média ponderada
sDM=DM[novoind]
berrDM = errDM[novoind] 
mediaDM[i] = weighted.mean(sDM,1/berrDM^2)
}
# fim das simulações

#estatísticas úteis:
summary(mediaDM)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.45   18.50   18.51   18.51   18.52   18.55
# quantis 0.1 e 0.9 dessas médias
# quantis são frequentemente usados para representar as incertezas nas medidas ou estimativas.
quantile(mediaDM, c(0.1, 0.9))
##      10%      90% 
## 18.48367 18.53086
# visualização dos resultados da simulação:
hist(mediaDM,col='salmon',xlab='DM do bootstrap',ylab='número',main='')

ATENÇÃO: deve-se sempre usar a média obtida diretamente dos dados, não do boostrap. Ele é útil para estimar erros e intervalos de confiança.

2. propagação de erros: relação entre fluxo e magnitude

Vamos fazer uma simulação para testar a fórmula da propagação de erros.

Vamos supor que uma fonte tem fluxo \(f = 10 \pm 1\) (em unidades arbitrárias)

Nesse caso diz-se que a razão sinal ruído é: 10/1 = 10

Vamos modelar a distribuição de fluxo como \(N(\mu = 10, \sigma = 1)\)

Podemos exprimir o fluxo em magnitudes como \(m=-2.5 \log f\) (ponto zero arbitrário).

Queremos estimar o erro das magnitudes.

# simulação dos fluxos
n=1000
# crio um vetor numérico para armazenar os valores de fluxo
x = numeric(n)
for(i in 1:n){
# como o fluxo deve ser positivo:
  x[i] = -1
  while(x[i] <= 0){x[i]=rnorm(1,mean=10,sd=1)}
}

mag=-2.5*log10(x)

hist(mag,col='cadetblue1',main='',ylab='número',breaks=50)

# desvio padrão das magnitudes
# (round arredonda o resultado, no caso com 3 casas decimais)
round(sd(mag),3)
## [1] 0.112

Compare o valor obtido da simuação com o valor esperado por propagação de erros: \[\sigma_m \simeq 1.0857 \frac{\sigma_f}{f}=0.109\]

3. otimização de funções uni-dimensionais

No caso geral, a minimização do \(\chi^2\) deve ser feita numericamente. O problema de se achar extremos de funções é denominado “otimização”. Há muitos tipos diferentes de problemas de otimização, como os sujeitos a vínculos, que não considerarei aqui.

Os problemas de maximização e minimização de funções são equivalentes. Vamos considerar apenas a minimização.

Vamos exemplificar com a função \(f(x) = x^{0.1}e^{-\sin [10 x \sin(x)]}\) entre 0 e 2.

# definição da função
f = function(x) {x^0.1*exp(-sin(10*x*sin(x)))}

# vetor com valores de x
x = seq(0,2,0.001)

# cálculo da variável dependente
y = f(x)

# visualização
plot(x,y, col='red',type='l',lwd=3)

Esta função tem vários mínimos neste intervalo. Vamos determinar a posição do mínimo próximo de 0.4.

optimize

A função optimize aplica-se apenas ao caso unidimensional e usa uma combinação de golden search e interpolações parabólicas. É preciso também dar dois valores para isolar o mínimo desejado (a seguir, 0.2 e 0.6). No exemplo abaixo vamos imprimir os valores de \(f(x)\) encontrados no processo de minimização:

f = function(x) {print(x^0.1*exp(-sin(10*x*sin(x))))} 
xmin = optimize(f, interval = c(0.2, 0.6))
## [1] 0.3524244
## [1] 0.3623167
## [1] 0.4162573
## [1] 0.3357775
## [1] 0.3357025
## [1] 0.3356316
## [1] 0.3356305
## [1] 0.3356305
## [1] 0.3356305
## [1] 0.3356305
xmin 
## $minimum
## [1] 0.3973244
## 
## $objective
## [1] 0.3356305

$minimum é a posição do mínimo, \(x_{min}\), e $objective é o valor de \(f\) neste mínimo, \(f(x_{min})\).

optim

A função optim (não confundir com optimize) serve tanto para otimizações em 1 dimensão como em \(N\) dimensões.

Para minimizações unidimensionais recomenda-se o método de Brent, que não envolve o cálculo de gradientes.

É preciso também dar um valor inicial (\(x_0\)) e dois valores para isolar o mínimo desejado:

# vamos reescrever a função de modo que ela não imprima os valores calculados em cada iteração:
f = function(x) {x^0.1*exp(-sin(10*x*sin(x)))} 
x0=0.6
optim(x0,f,method='Brent',lower=0.2,upper=0.6)
## $par
## [1] 0.397314
## 
## $value
## [1] 0.3356305
## 
## $counts
## function gradient 
##       NA       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Nesse caso a função retorna o valor do mínimo e o da função nesse ponto.

$convergence = 0 significa que a minimização convergiu.

Uma opção envolvendo o cálculo de derivadas é o CG: gradientes conjugados. Nesse caso só o valor inicial da busca deve ser fornecido e o algoritmo retorna o mínimo mais próximo a esse valor:

optim(x0,f,method='CG')
## $par
## [1] 0.3973139
## 
## $value
## [1] 0.3356305
## 
## $counts
## function gradient 
##      135       18 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Para se obter apenas o valor do mínimo:

xmin = optim(0.4,f,method='CG')$par
xmin
## [1] 0.3973135

Há inúmeros algoritmos de minimização e alguns podem ser melhores que outros em certos casos. Faça ?optim e ?optimize para conhecer melhor essas funções.

4. otimização de funções bi-dimensionais

Consideremos, inicialmente, a função de Rosenbrok (função banana): \[ f(x,y) = (1 − x)^2 + 100 (y − x^2)^2 \] que tem um mínimo em (1,1).

f = function(x,y) {return((1-x)^2+100*(y-x^2)^2)}
x = seq(-3,3,0.2)
y = seq(-3,3,0.2)
z = outer(x,y,f)
persp(x,y,z,phi=45,theta=-10,col="yellow",shade=.00000001,ticktype="detailed", main = 'função banana')

fl = function(x,y) {return(log((1-x)^2+100*(y-x^2)^2))}
z = outer(x,y,fl)
image(x,y,z, col = rainbow(100), main = 'log da função banana')

otimização com optim

Vamos reescrever \(f\) para que os parâmetros entrem como um vetor, \(x\). Vamos considerar o ponto (0,0) como ponto inicial.

Note que para dimensões maiores que 1 não se isola o mínimo mas define-se um ponto inicial para começar a minimização (que pode convergir para um mínimo secundário!).

f = function(x) {(1-x[1])^2 + 100*(x[2]-x[1]^2)^2}
# valor inicial dos parâmetros
xini=c(0,0)
# valor estimado
optim(xini,f, method = "Nelder-Mead")$par
## [1] 0.9999564 0.9999085

A função de Himmelblau tem 4 mínimos em (−3.7793,−3.2832); (−2.8051, 3.1313); (3,2); (3.5844, −1.8481 ): \[ f (x, y) = (x^2 + y − 11)^2 + (x + y^2 − 7)^2 \]

f = function(x1,y1) {(x1^2 + y1 - 11)^2 + (x1 + y1^2 - 7)^2}
x <- seq(-4.5,4.5,by=.4)
y <- seq(-4.5,4.5,by=.4)
z <- outer(x,y,f)
persp(x,y,z,phi=45,theta=45,col="yellow",shade=.65,ticktype="detailed")

image(x,y,z, col = rainbow(100))
# mínimos
points(x=-3.7793,y=-3.2832,pch=20)
points(x=-2.8051,y= 3.1313,pch=20)
points(x=3.5844,y= -1.8481,pch=20)
points(x=3,y=2,pch=20)

Vamos variar os métodos e o ponto inicial:

f = function(x) (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
# escrevendo a função com $par no final retorna apenas o valos dos parâmetros (x[1],x[2]) no mínimo
optim(c(-4,-4), f, method = 'Nelder-Mead')$par
## [1] -3.779347 -3.283172
optim(c(-3.8,-4.1), f, method = 'BFGS')$par
## [1] -3.779310 -3.283186
optim(c(2,-2), f, method = 'CG')$par     
## [1]  3.584428 -1.848126
optim(c(-4,4), f, method = 'L-BFGS-B')$par  
## [1] -2.805118  3.131312
optim(c(2,2), f, method = 'L-BFGS-B')$par  
## [1] 3 2
optim(c(2,2), f, method = 'Nelder-Mead')$par
## [1] 3.000014 2.000032

Opções da função optim:

# note que o method é opcional; o default é Nelder-Mead (simplex)
r = optim(c(2,2), f)
#checar a convergência
r$convergence == 0 # TRUE se convergiu
## [1] TRUE
#solução
r$par
## [1] 3.000014 2.000032
#valor da função no mínimo
r$value
## [1] 3.37533e-08

Mas nem sempre a vida é fácil: como a função optim pode ser usada em casos multidimensionais, vamos considerar uma hiper-banana em 6 dimensões que tem um mínimo em (1,1,1,1,1,1):

f = function(x) {(1-x[1])^2 + 100*(x[2]-x[1]^2)^2+ 33*(x[2]-x[3])^2+120*(1-x[4])^2+44*(1-x[5]^2)^2+ 133*(x[6]-x[1])^2}
xini=c(0,0,0,0,0,0)
res = optim(xini,f, method = "Nelder-Mead")

#convergência: o valor de 'convergence' deve ser 0
res$convergence
## [1] 1
res$par
## [1] 0.25480283 0.04798550 0.04785831 1.01955751 1.00382755 0.24865746
res$value
## [1] 0.637523
#como não convergiu, vamos mudar o metodo:

res = optim(xini,f, method = "L-BFGS-B")
res$convergence
## [1] 0
res$par
## [1] 0.9998721 0.9997361 0.9997219 0.9999996 0.0000000 0.9998768
res$value
## [1] 44
#convergiu para um mínimo errado! vamos recomeçar com o último ponto e com outro método:
xini=res$par
res = optim(xini,f, method = "Nelder-Mead")
res
## $par
## [1] 0.9492360 0.9006119 0.8997636 0.9999170 1.0000543 0.9495540
## 
## $value
## [1] 0.002634632
## 
## $counts
## function gradient 
##      435       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Ufa, agora deu certo!!

Moral da história: problemas diferentes podem precisar de métodos diferentes!

Em muitos casos é necessário se repetir o processo com diferentes inicializações!

Estes métodos podem ser usados em otimizações de funções com mais de duas dimensões.

5. regressão linear

Vamos gerar uma amostra \((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:

# reprodutibilidade
set.seed(12345)

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

estimativa dos parâmetros de um modelo linear

Vamos implementar as relações para a obtenção dos parâmetros de uma reta: \[ 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} \]

# programação convencional (tipo fortran, c, ...) usando loops
S=0
Sx=0
Sy=0
Sxx=0
Sxy=0
for (i in 1:n) {
S=S+1./err[i]^2
Sx=Sx+x[i]/err[i]^2
Sy=Sy+y[i]/err[i]^2
Sxx=Sxx+x[i]^2/err[i]^2
Sxy=Sxy+x[i]*y[i]/err[i]^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] 0.8560904
## [1] 1.992887
#erros nos parâmetros
siga; sigb
## [1] 0.2709611
## [1] 0.03731369
#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)

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

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

legend('topleft',c('simulação','ajuste'),col=c('green','black'),pch = c(20,20))

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

# chi2
chi2=0
for (i in 1:n) {
  chi2=chi2+((y[i]-a-b*x[i])/err[i])^2
}
chi2
## [1] 8.737462
# chi2 reduzido
nu=n-2
chi2/nu
## [1] 1.092183

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.

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!

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  
##      0.8561       1.9929

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

Comparem os coeficientes calculados com os usados na simulação dos dados: a diferença é devida ao ruído.

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  
##      0.8561       1.9929
#summary
summary(ajuste)
## 
## Call:
## lm(formula = y ~ x, weights = 1/err^2)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9325 -0.7483 -0.3806  1.1706  1.2972 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.8561     0.2832   3.023   0.0165 *  
## x             1.9929     0.0390  51.105 2.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.045 on 8 degrees of freedom
## Multiple R-squared:  0.9969, Adjusted R-squared:  0.9966 
## F-statistic:  2612 on 1 and 8 DF,  p-value: 2.381e-11
#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 a relação é significativa

# 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) 
##   0.8560904
##        x 
## 1.992887
# essas quantidades podem ser obtidas também com algumas funções:
# coeficientes ajustados:
coefficients(ajuste)
## (Intercept)           x 
##   0.8560904   1.9928869
# resíduos do ajuste:
residuals(ajuste)
##           1           2           3           4           5           6 
## -3.05149352  1.37152858 -0.27202795 -0.27303190 -1.46216529 -0.04758527 
##           7           8           9          10 
##  3.44975766  0.86537199  1.15794253 -1.10144691
# valores ajustados para cada x:
fitted.values(ajuste)
##         1         2         3         4         5         6         7         8 
## 22.406290 27.035845 23.604366 27.345281 14.501814  5.829493 10.574266 16.078488 
##         9        10 
## 22.609605 30.442597

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

modelos em lei de potência

Modelos em lei de potência podem ser tratados como lineares. Por exemplo, vamos supor que \(y = ax^{b}\) para \(a > 0\) e \(x > 0\). Logo, \(\log y = \log a + b \log x\), que tem a forma conhecida de uma reta. Vamos ilustrar este tipo de modelo com a função de massa inicial.

A Função de Massa Inicial (IMF) dá a distribuição em massa das estrelas formadas em um surto de formação estelar:

\(\xi (m) dm\): número de estrelas formadas com massa entre \(m\) e \(m+dm\)

Os dados abaixo são da Tabela 2 do artigo de Salpeter (1955, ApJ, 121, 161), onde ele propõe um modelo que até hoje é conhecido como modelo de Salpeter para a IMF:

# x = log(m/M_sun) + 1
x = c(2.23,2.08,1.93,1.78,1.63,1.48,1.33,1.20,1.09,1.,0.93,0.86,0.80,0.74,0.68,0.62,0.56,0.50)

# x = log(m/M_sun)
x = x-1

# y = log xi + 10
y = c(6.63,7.10,7.36,7.52,7.72,8.00,7.98,7.98,8.32,8.50,8.60,8.70,8.83,8.97,9.04,9.13,9.20,9.22)

# número de dados
length(x); length(y)
## [1] 18
## [1] 18
# modelo linear considerando todos os pontos
modelo=lm(y~x)

plot(x,y,main='IMF de Salpeter', xlab='log(m/M_sun)', ylab='log xi + 10',pch=20,cex=2,col='red')
abline(modelo,col='blue',lwd=2)

# modelo ajustado:
modelo
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##       8.535       -1.404
# no artigo, Salpeter considera apenas os pontos em -0.4 < log(m/M_sun) < 1 para fazer o ajuste:
xl = x[x > -0.4 & x < 1]
yl = y[x > -0.4 & x < 1]

# número de dados neste ajuste
length(xl)
## [1] 14
# modelo linear considerando apenas os pontos em -0.4 < log(m/M_sun) < 1 
modelo=lm(yl~xl)
plot(x,y,main='IMF de Salpeter', xlab='log(m/M_sun)', ylab='log xi + 10',pch=20,cex=2,col='red')
points(xl,yl,col='blue',cex=2,pch=20)
abline(modelo,col='red',lwd=2)

# modelo ajustado:
modelo
## 
## Call:
## lm(formula = yl ~ xl)
## 
## Coefficients:
## (Intercept)           xl  
##       8.531       -1.346

Com esse resultado ele propõe o que hoje chamamos de função de massa inicial de Salpeter: \[ \xi (m) \propto m^{-1.35} \]

6. mínimos quadrados linear generalizado

Consideremos uma função \[ y(x) = a_1 + a_2 X_2(x) + a_3 X_3(x)+ ...\] onde os \(X_k(x)\) são funções arbitrárias de \(x\). Note que a função é linear com respeito aos parâmetros \(a_k\).

Para um certo conjunto de dados \(\{x_i,y_i,\sigma_i \}\), os parâmetros \(a_k\) podem ser obtidos minimizando o \(\chi^2\): \[ \chi^2 = \sum_{i=1}^{N} \Bigg[ {y_i - \sum_{k=1}^{M} a_k X_k(x_i) \over \sigma_i} \Bigg]^2 \]

ou, em linguagem matricial, \[ \chi^2 = \mathbf{(b-A.a)^T . (b-A.a)} \]

onde \[ A_{ij} = {X_j(x_i) \over \sigma_i}~~~~~~~~~b_i = {y_i \over \sigma_i}\]

a solução para os coeficientes é: \[ a_j = \sum_{k=1}^{m} [\alpha]_{jk}^{-1} \beta_k = \sum_{k=1}^{m} C_{jk} \beta_k,\] onde
\[ \alpha_{kj} = \sum_{i=1}^n {X_j(x_i) X_k(x_i) \over \sigma_i^2}~~~~~~~~~~ \beta_k = \sum_{i=1}^n {y_i X_k(x_i) \over \sigma_i^2} \]
C é a matriz de covariância, \(C = \alpha^{-1}\)

os erros dos parâmetros são dados por \(\sigma_j^2 = C_{jj}\)

Vamos considerar N dados gerados de acordo com a função \[ y= a_1 + a_2 x + a_3 x^2 + N(0,sd=\sigma) \] isto é, um polinômio de segundo grau mais ruído gaussiano:

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

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

Para obter os coeficientes do ajuste comecemos por construir as matrizes relevantes (vamos usar uma programação convencional):

X=matrix(rep(0,N*M),N)
for (i in 1:N) {
for (j in 1:M) {
X[i,j]= x[i]^(j-1)
}
}

alpha=matrix(rep(0,M*M),M)
for (j in 1:M) {
for (k in 1:M) {
for (i in 1:N) {
alpha[j,k]=alpha[j,k]+X[i,j]*X[i,k]/sig[i]^2
}
}
}

beta=rep(0,M)
for (k in 1:M) {
for (i in 1:N) {
beta[k]=beta[k]+y[i]*X[i,k]/sig[i]^2
}
}

A matriz de covariância é a inversa de alpha. Antes é bom ver se o determinante de alpha é diferente de zero (se não tem que usar SVD: singular value decomposition)

det(alpha)
## [1] 4034966

Matriz inversa:

C=solve(alpha)

teste de sanidade- a matriz multiplicada por sua inversa deve ser igual à matriz unidade (note o símbolo para multiplicação matricial):

U=C%*%alpha
U
##               [,1]         [,2]          [,3]
## [1,]  1.000000e+00 1.776357e-14  1.421085e-14
## [2,] -3.552714e-15 1.000000e+00 -2.842171e-14
## [3,]  4.440892e-16 1.776357e-15  1.000000e+00

Os parâmetros então podem ser calculados como:

a=rep(0,M)
for (j in 1:M) {
for (k in 1:M) {
a[j]=a[j]+C[j,k]*beta[k]
}
}
yp=a[1]+a[2]*xp+a[3]*xp^2
plot(x,y,pch=20,cex=2,main = 'v: modelo   a: ajuste')
lines(xp,yp,col='blue',lwd=3)
lines(xm,ym,col='red',lwd=3)

# parâmetros:
a
## [1]  0.5719313 -0.4525223  0.4220399

erros nos parâmetros:

erros=rep(0,M)
for (j in 1:M) { 
  erros[j] = sqrt(C[j,j])
}
erros
## [1] 0.26986020 0.29773624 0.06232058

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

Vamos agora aplicar alguns dos métodos das seções anteriores.

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.

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)

# 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

relação de Faber-Jackson: relação entre a luminosidade (ou magnitude absoluta) e a dispersão de velocidades

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

plot(dados$log_sigma,dados$BT,xlab=expression(paste('log ',sigma,' (km/s)')),ylab='BT mag',main='aglomerado de Virgo',pch = 20, cex=2)

# ajuste com lm()
fj = lm(dados$BT ~ dados$log_sigma)
abline(fj,col='red',lwd=3)

summary(fj)
## 
## Call:
## lm(formula = dados$BT ~ dados$log_sigma)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05711 -0.34011  0.01854  0.28725  1.42640 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      23.0618     1.3464  17.128 1.37e-12 ***
## dados$log_sigma  -5.2276     0.6047  -8.645 7.98e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.628 on 18 degrees of freedom
## Multiple R-squared:  0.8059, Adjusted R-squared:  0.7951 
## F-statistic: 74.74 on 1 and 18 DF,  p-value: 7.98e-08

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)

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

plot(dados$log_Ae,dados$Sigmae,xlab='log Ae (arcsec)',ylab=expression(paste(Sigma,' (mag/arcsec^2)')),main='aglomerado de Virgo',pch = 20, cex=2)

# ajuste
rk = lm(dados$Sigmae ~ dados$log_Ae)
abline(rk,col='red',lwd=3)

summary(rk)
## 
## Call:
## lm(formula = dados$Sigmae ~ dados$log_Ae)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7809 -0.5433 -0.1028  0.4030  1.2337 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   18.2793     0.7783  23.486 5.91e-15 ***
## dados$log_Ae   1.5083     0.4234   3.562  0.00223 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6521 on 18 degrees of freedom
## Multiple R-squared:  0.4135, Adjusted R-squared:  0.3809 
## F-statistic: 12.69 on 1 and 18 DF,  p-value: 0.002228

Dá para notar que esta relação é muito mais ruidosa que a de Faber-Jackson.

Vamos agora considerar o “Plano Fundamental”: um melhor ajuste é obtido combinando-se 3 variáveis: \(BT, \sigma, \Sigma_e\)

Vamos considerar \(BT = f(\Sigma_e,\sigma)\):

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

# a função lm() pode ser usada para 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)

library(FactoClass)
## Loading required package: ade4
## Loading required package: ggplot2
## Loading required package: ggrepel
## Loading required package: xtable

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\).

relação Dn-sigma- um outro indicador de distância

Dn: diâmetro da abertura circular (em arcsec) contendo o brilho superficial integrado de 20.75 mag/arcsec^2

par(mfrow=c(1,1))
plot(dados$log_sigma,dados$log_Dn,xlab=expression(paste('log ',sigma,' (km/s)')),ylab='log Dn (arcsec)',main='aglomerado de Virgo',pch=20,cex=2)
ds = lm(dados$log_Dn ~ dados$log_sigma)
abline(ds,col='red',lwd=2)

summary(ds)
## 
## Call:
## lm(formula = dados$log_Dn ~ dados$log_sigma)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.156787 -0.048838  0.001943  0.049526  0.111214 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.75978    0.15871  -4.787 0.000147 ***
## dados$log_sigma  1.12293    0.07128  15.755 5.65e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07402 on 18 degrees of freedom
## Multiple R-squared:  0.9324, Adjusted R-squared:  0.9286 
## F-statistic: 248.2 on 1 and 18 DF,  p-value: 5.654e-12

comparação do resultado dos 4 ajustes com o \(\chi_{red}^2\):

sum(fj$residuals^2)/(fj$df.residual)
## [1] 0.394336
sum(rk$residuals^2)/(rk$df.residual)
## [1] 0.4252518
sum(pf$residuals^2)/(pf$df.residual)
## [1] 0.1036415
sum(ds$residuals^2)/(ds$df.residual)
## [1] 0.005479083

vemos que o Plano Fundamental é um indicador de distância melhor que as relações de Faber-Jackson e Kormendy, e que Dn-sigma tem um desempenho melhor que o do Plano Fundamental.

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

Vamos agora voltar ao problema específico da minimização do \(\chi^2(\mathbf{a})\) para obtenção dos parâmetros \(\mathbf{a}\) e de seus erros: \[ \chi^2 (\mathbf{a}) = \sum_{i=1}^N {(y_i - f(x_i;\mathbf{a}))^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] \] com 3 parâmetros: \(\mathbf{a} = \{\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:
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)

# 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.64262    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: 3.619e-07
# 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

# transformando a tabela em raio e brilho superficial
raio=radius
bs=surf_mag

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

# função chi2 
fchi2 = function(x) {
s=0
for (k in 1:npt) {
r=raio[k]
bso=bs[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, 
   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)*((raio/res$par[2])^{1/res$par[3]}-1)
lines(raio,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:

# transformando a tabela em raio e brilho superficial
raio=radius
bs=surf_mag

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

# 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=raio[d[k]]
bso=bs[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.08355144
sd(p2)
## [1] 10.23274
sd(p3)
## [1] 0.1262216

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:   
cor.test(p1,p2,method='spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  p1 and p2
## S = 1282808, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9923031
cor.test(p1,p3,method='spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  p1 and p3
## S = 7224722, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9566516
cor.test(p2,p3,method='spearman')
## 
##  Spearman's rank correlation rho
## 
## data:  p2 and p3
## S = 12104458, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9273732