aula de hoje:

    1. combinação de medidas
    1. propagação dos erros
    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. exercícios
  • Apêndice: otimização de funções

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 usar bootstrap para estimar uma barra de erro para esta estimativa: vamos reamostrar o vetor DM muitas vezes, com substituição, calculando a média ponderada cada vez, e determinar o desvio padrão e outras estatísticas da distribuição das médias ponderadas simuladas:

# 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]
serrDM = errDM[novoind] 
mediaDM[i] = weighted.mean(sDM,1/serrDM^2)
}
# fim das simulações

# visualização dos resultados da simulação:
hist(mediaDM,col='salmon',xlab='DM do bootstrap',ylab='número',main='')
# valor médio dos dados:
abline(v=DM_MV,lwd=3)

#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
# desvio padrão:
sd(mediaDM)
## [1] 0.01857035
# opção robusta para o desvio padrão:
sig_G = 0.7413*(quantile(mediaDM,0.75,names = FALSE) - quantile(mediaDM,0.25,names = FALSE))
sig_G
## [1] 0.01857339
# no caso o desvio padrão e sigma_G são muito parecidos!

# intervalo de confiança de 95%:
quantile(mediaDM, c(0.025, 0.975))
##     2.5%    97.5% 
## 18.47067 18.54440

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

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

Fluxos e magnitudes se relacionam como: \[m=m_0-2.5 \log f\]

Se o fluxo medido tem um erro \(\sigma_f\), qual será o erro na magnitude?

Vamos usar a fórmula da propagação de erros. No caso, \[\sigma_m^2 = \Big({dm \over df}\Big)^2 \sigma_f^2\] e \[{dm \over df} = \Bigg|{-2.5/ln10 \over f}\Bigg| \simeq {1.086 \over f}\]

Logo, \[\sigma_m \simeq 1.086 {\sigma_f \over f}\]

Vamos fazer uma simulação para testar este resultado.

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

Vamos supor aqui que \(m_0=0\) e, portanto, \(m=-2.5 \log f\).

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

# magnitudes
mag=-2.5*log10(x)

# distribuição das magnitudes
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 do desvio padrão obtido da simulação com o valor esperado por propagação de erros: \[\sigma_m \simeq 1.0857 \frac{\sigma_f}{f}=0.109\]

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

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} \]

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

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 (exatamente o que fizemos mais acima):

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)

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 é uma reta na variável \(\log x\). 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)

# 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
# valor do expoente:
modelo$coefficients[2]
##        xl 
## -1.345558

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

Por esse modelo, em um surto de formação estelar a razão entre o número de estrelas de \(10 M_\odot\) em relação às de \(1 M_\odot\) é \(10^{-1.35} \simeq 0.044\)

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

4. 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 (\{a_k\}) = \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 (\{a_k\}) = \mathbf{(b-A.a)^T . (b-A.a)} \]

onde \(\mathbf{a} = \{a_1, a_2, ...\}\) e \[ 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}\)

Exemplo: 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 é, o polinômio de segundo grau mais ruído gaussiano discutido anteriormente.

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] 4908316

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 3.552714e-14 1.136868e-13
## [2,]  0.000000e+00 1.000000e+00 0.000000e+00
## [3,] -4.440892e-16 0.000000e+00 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]*xm+a[3]*xm^2
plot(x,y,pch=20,cex=2,main = 'modelo e ajuste')
lines(xm,yp,col='blue',lwd=3)
lines(xm,ym,col='red',lwd=3)
legend('topleft',c('modelo','ajuste'),col=c('red','blue'),pch = c(20,20))

# parâmetros:
a
## [1]  1.3721938 -1.2047424  0.5138822

erros nos parâmetros:

erros=rep(0,M)
for (j in 1:M) { 
  erros[j] = sqrt(C[j,j])
}
erros
## [1] 0.34337813 0.30014698 0.05516364

5. 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 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\]

vamos usar lm() para determinar esta relação:

plot(dados$log_sigma,dados$BT,xlab=expression(paste('log ',sigma,' (km/s)')),ylab='BT mag',main='Virgo: relação de Faber-Jackson',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

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!

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

# ajuste com lm()
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)\) 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\).

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='Virgo: relação Dn-sigma',pch=20,cex=2)

# ajuste com lm()
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.

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

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

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

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

7. exercícios

  1. inicialize a semente de números aleatórios com seu número USP.
  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. Discuta se as propriedades químicas das galáxias early-type, descritas pelo Mg2, podem ser estimadas do brilho superficial e da dispersão de velocidades. Use mínimos quadrados para analisar a relação do Mg2 com o brilho superficial e a dispersão de velocidades, tanto com cada variável separadamente quanto com as duas juntas. Escreva suas conclusões.
  1. Use bootstrap para estimar o erro no expoente da função de massa de Salpeter. Como essas incertezas afetam a razão entre o número de estrelas formadas com 10 \(M_\odot\) em relação às com 1 \(M_\odot\)? Discuta a distribuição dessa razão com os resultados do bootstrap.

Apêndice: otimização de funções

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

Exercício: encontre o mínimo próximo a 1 da função acima.

A2. 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(1000), main = 'log da função banana')

# agora usando a biblioteca fields
suppressMessages(library(fields))
z = outer(x,y,f)
image.plot(x, y, z)

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
# compare com o esperado: (1,1)!

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,cex=2)
points(x=-2.8051,y= 3.1313,pch=20,cex=2)
points(x=3.5844,y= -1.8481,pch=20,cex=2)
points(x=3,y=2,pch=20,cex=2)

# ou
image.plot(x, y, z)
points(x=-3.7793,y=-3.2832,pch=20,cex=2,col='salmon')
points(x=-2.8051,y= 3.1313,pch=20,cex=2,col='salmon')
points(x=3.5844,y= -1.8481,pch=20,cex=2,col='salmon')
points(x=3,y=2,pch=20,cex=2,col='salmon')

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)
r
## $par
## [1] 3.000014 2.000032
## 
## $value
## [1] 3.37533e-08
## 
## $counts
## function gradient 
##       59       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
#cheque a convergência: $convergence == 0 - convergiu
#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.