Continuando o assunto de probabilidade, vamos ver o problema de mínimos quadrados linear como um problema de maximização da verossimilhança.
Os dois métodos são equivalentes, quando a distribuição de probabilidade do erro das amostras é normal.
Mas a verossimilhança envolve uma interpretação diferente do problema, com aspectos probabilísticos.
Como consequência, podemos quantificar as incertezas na escolha dos parâmetros e nas predições do modelo.
using Distributions
using Plots
using Random
using LinearAlgebra: ⋅
De um ponto de vista probabilístico, consideramos incertezas inerentes na obtenção dos dados $x$ e na definição do modelo.
Com base nisso, interpretamos o modelo como uma relação para o valor esperado de $y$, dado $x$, i.e.
Além disso, assumimos que os erros são independentes entre si. Ou seja, em quaisquer dois pontos $x_i, x_j$, os erros $\epsilon_i, \epsilon_j$ nesses pontos são independentes entre si.
Como os erros são normais, serem independentes é equivalente a não estarem correlacionados, o que pode ser expresso por $E(\epsilon_i\epsilon_j)=0$.
Observe que quanto mais perto $\beta_0 + \beta_1 x$ estiver de $y$, maior será a verossimilhança.
Mas independentemente da forma da função de densidade de distribuição, maximizar a verossimilhança é aumentar a probabilidade do modelo resultar no dado observado.
Dado um conjunto de observações (independentes) $(\mathbf{x}, \mathbf{y}) = (x_i, y_i)_{i=1}^N$, a probabilidade conjunta é o produto das probabilidades, $\mathcal{P}(y_1|x_1)\mathcal{P}(y_2|x_2)\cdots\mathcal{P}(y_N|x_N)$.
Em relação à função densidade de probabilidades, obtemos
Como o logaritmo é monónoto crescente, maximizar a verossimilhança é equivalente a maximizar a log-verossimilhança.
No caso da normal, vemos que a log-verossimilhança tem uma dependência bem mais amigável, no sentido de ser linear.
No caso de um conjunto de dados, a função de log-verossimilhança toma a forma
O primeiro termo não depende de $\boldsymbol\beta$.
E vemos, do segundo termo, que maximizar a verossimilhança é equivalente a minimizar o erro quadrático (i.e. a soma dos quadrados dos resíduos (RSS))
Com esse arcabouço probabilístico, podemos extrair informações sobre a incerteza na determinação dos parâmetros $\boldsymbol\beta=(\beta_0, \beta_1)$.
Obtivemos, acima, o maximizador da verossimilhança $\hat{\boldsymbol\beta}$ dado exatamente pela forma normal da solução pelo método de mínimos quadrados:
Nesta fórmula, $X = [\mathbf{1}, \mathbf{x}]$ é a matrix de Vandermonde, que assumimos de posto máximo, e $\mathbf{y}$ é dado por $\mathbf{y} = (y_1, \ldots, y_N)$
Agora, ao invés de olharmos estritamente para $\boldsymbol\beta$ que minimiza o erro entre as medições $y_i$ e o do resultado do modelo $\beta_0 + \beta_1 x_i$, vamos olhar para possíveis outras escolhas $\mathbf{y} = X\mathbf{\boldsymbol\beta} + \boldsymbol{\epsilon}$, que talvez extrapolassem melhor $y = \beta_0 + \beta_1 x$ para outros dados.
Considerando, então, $\mathbf{y} = X\mathbf{\boldsymbol\beta} + \boldsymbol{\epsilon}$, obtemos
Se $X$ e $\epsilon$ fossem escalares, seria fácil deduzir que $\operatorname{Var}(\beta) = (X^TX)^{-1}X^T\operatorname{Var}(\epsilon)X(X^TX)^{-1}$. Mas não é o caso.
No caso multidimensional, em que $X$ é de fato uma matriz e $\boldsymbol\epsilon$ e $\boldsymbol\beta$ são vetores, temos o valor esperado nos dando um vetor com os valores esperados de cada coordenada e temos a variância sendo estendida a uma matrix de variância-covariância, nos dando não apenas a variação de cada coordenada, mas também uma variação conjunta entre cada par de coordenadas. E a variância de $\boldsymbol\beta$ envolve tudo isso.
Naturalmente, $\operatorname{Cov}(\beta_i, \beta_i) = \operatorname{Var}(\beta_i)$.
Juntando os dois conceitos, temos a matriz de variância-covariância
Estamos interessados, então, na matriz $\operatorname{Cov}(\boldsymbol\beta)$, sabendo que $\boldsymbol\beta - \hat{\boldsymbol\beta} = - (X^TX)^{-1}X^T\boldsymbol{\epsilon}$.
Nesse caso, usamos que
onde $I$ é a matriz identidade.
Com a incerteza nos parâmetros caracterizada pela matriz de variância-covariância $\operatorname{Cov}(\boldsymbol\beta)$, podemos estimar a incerteza do modelo
Considerando $y = \beta_0 + \beta_1 x$, temos
A questão é como calcular essa variância do lado direito da expressão.
Explicitando a expressão e usando que $E(\beta_0 + \beta_1 x) = \hat\beta_0 + \hat\beta_1 x$, temos
O processo de determinação dos parâmetros envolve uma média que faz com que variância calculada acima seja, de fato, uma medida da incerteza no parâmetro.
De fato, observemos o caso mais simples em que queremos modelar apenas os erros $\epsilon = N(0,\sigma)$ por uma constante.
Fazemos uma série de amostras $(\epsilon_i)_{i=1}^N$ e buscamos encontrar o valor médio $\hat\beta_0$ pelo método de mínimos quadrados.
Nesse caso, a matrix $X$ é uma matrix $N\times 1$ com todos os elementos iguais a 1:
Observe que essas medidas dependem do desvio padrão $\sigma$ da incerteza na coleta dos dados.
Em certos casos, essa medida pode ser obtida dos instrumentos de medição e do próprio processo de coleta.
Na falta de maiores informações, uma alternativa é utilizar o desvio padrão corrigido da amostra. No caso, no entanto, não estamos fazendo um conjunto $(x_j,y_j)$ de medições na mesma situação. Estamos medindo em pontos diversos.
Por conta disso, usamos o próprio modelo para compensar isso e devemos considerar que agora o grau de liberdade disponível é $N-m$, onde $m=2$.
No caso de dados nas "mesmas" condições, temos $m=1$, como feito antes. Mas agora, temos dois parâmetros, $\beta_0$ e $\beta_1$, logo $m=2$. Em outros modelos, $m$ pode ser maior.
Assim, o desvio padrão corrigido da amostra é
Vamos construir um exemplo sintético.
Vamos "perturbar" uma determinada reta com erros aleatórios segundo uma determinada normal com desvio padrão pré-definido.
b = 1.0
m = 0.2
σ = 1.0
N = 20
20
Random.seed!(1200)
x_org = [0.0, 10.0]
y_org = b .+ m * x_org
data_x = collect(range(0.0,10.0, length=N)) + rand(N)/N
data_y = b .+ m * data_x .+ rand(Normal(0,σ), N)
plot(x_org, y_org, label = "y = b + mx", title = "Dados perturbados aleatoriamente, em um exemplo sintético", titlefont = 10)
scatter!(data_x, data_y, markersize=3, label="amostra")
X = [ones(N) data_x]
β = X \ data_y
2-element Vector{Float64}: 1.1988117103655653 0.1431233525841148
plot(x_org, y_org, label="original", legend=:topleft, title = "Ajuste aos dados", titlefont = 10)
scatter!(data_x, data_y, markersize=3, label="amostra")
plot!(x_org, β[1] .+ β[2] * x_org, label="modelo")
Cov_β = σ^2 * inv(X' * X)
2×2 Matrix{Float64}: 0.187317 -0.0273258 -0.0273258 0.0054378
println("Valores originais: β₀=$b, β₁=$m")
println("CI 68% para β₀: [$(round(β[1] - Cov_β[1,1],digits=2)), $(round(β[1] + Cov_β[1,1],digits=2))]")
println("CI 68% para β₁: [$(round(β[2] - Cov_β[2,2],digits=2)), $(round(β[2] + Cov_β[2,2],digits=2))]")
println("CI 95% para β₀: [$(round(β[1] - 2Cov_β[1,1],digits=2)), $(round(β[1] + 2Cov_β[1,1],digits=2))]")
println("CI 95% para β₁: [$(round(β[2] - 2Cov_β[2,2],digits=2)), $(round(β[2] + 2Cov_β[2,2],digits=2))]")
Valores originais: β₀=1.0, β₁=0.2 CI 68% para β₀: [1.01, 1.39] CI 68% para β₁: [0.14, 0.15] CI 95% para β₀: [0.82, 1.57] CI 95% para β₁: [0.13, 0.15]
Var_y = [[1; x] ⋅ (Cov_β * [1; x]) for x in data_x]
20-element Vector{Float64}: 0.18506906749986532 0.1598972882552096 0.13472605333885954 0.1138367091305883 0.09548818390233753 0.08014032619999441 0.0683046673239387 0.059364145036716144 0.053257487510014974 0.05035707669772974 0.05032645779142901 0.053546214207226044 0.059513423230789925 0.06796962535835986 0.08049659618144059 0.09529307440702718 0.11444701290418924 0.13374704790809286 0.15834618382131746 0.18587335929487245
plot(x_org, y_org, label="original", size=(600,300), titlefont=10,
title="Reta original, amostra, modelo, CI 68% (sombreado), CI 95% (barras de erro)", )
scatter!(data_x, data_y, markersize=3, label="amostra")
plot!(data_x,β[1] .+ β[2] * data_x, grid=false, yerror=2*sqrt.(Var_y),
ribbon=sqrt.(Var_y), fillalpha=0.2, label="modelo", legend=:topleft)
Como a amostra foi obtida de maneira sintética, conhecemos o desvio padrão associado a distribuição dos erros.
Mas em geral isso não está disponível.
Nesse caso, podemos aproximá-lo com o desvio padrão corrigido $s_q$ da amostra.
Observe que $s_q$ fica bem próximo do valor original de $\sigma$ e que o desvio padrão não-corrigido da amostra fica ligeiramente distante.
s_q = √(sum(abs2, data_y .- β[1] .- β[2] * data_x)/(N-2))
println("Desvio padrão original: $σ")
println("Desvio padrão corrigido da amostra: $(round(s_q, digits=3))")
println("Desvio padrão não corrigido da amostra: $(round(s_q * √((N-2)/N),digits=3))")
Desvio padrão original: 1.0 Desvio padrão corrigido da amostra: 0.953 Desvio padrão não corrigido da amostra: 0.904
Cov_q = s_q^2 * inv(X' * X)
Var_q = [[1; x] ⋅ (Cov_q * [1; x]) for x in data_x]
20-element Vector{Float64}: 0.16813890537256826 0.14526984645498642 0.12240128207042794 0.10342290001783515 0.0867529022231941 0.07280907017850197 0.06205614017962694 0.05393349898874427 0.04838547993890799 0.045750399394674124 0.045722581513173675 0.04864779384947513 0.054069121178268606 0.06175174793241178 0.07313274849766525 0.08657564139072518 0.10397737305768065 0.12151183628833889 0.14386063742207578 0.16886962036367678
plot(x_org, y_org, label="original", size=(600,300), titlefont=10,
title="Reta original, amostra, modelo, CI 68% (sombreado), CI 95% (barras de erro)\nusando o desvio padrão corrigido", )
scatter!(data_x, data_y, markersize=3, label="amostra")
plot!(data_x,β[1] .+ β[2] * data_x, grid=false, yerror=2*sqrt.(Var_q),
ribbon=sqrt.(Var_q), fillalpha=0.2, label="modelo", legend=:topleft)
Vamos, agora, refazer o exemplo da relação entre o salário anual médio e o grau de instrução, dos EUA.
Dados obtidos de US Bureau of Labor Statistics: Learn more, earn more: Education leads to higher wages, lower unemployment.
Nível de instrução | Média de salário semanal (USD) | Taxa de desemprego (%) |
---|---|---|
Doutorado | 1883 | 1,1 |
Profissional | 1861 | 1,6 |
Mestrado | 1497 | 2,0 |
Graduação | 1248 | 2,2 |
Associado* | 887 | 2,7 |
Graduação incompleta | 833 | 3.3 |
Ensino Médio | 746 | 3,7 |
Ensino Fundamental | 592 | 5,4 |
data_y = [592, 746, 833, 887, 1248, 1497, 1861, 1883]
plot(data_y, seriestype = :scatter, xlims=(0,9), ylims=(0,2000),
xticks=0:9, xaxis = "nível de instrução", yaxis="salário (USD)",
title="Média salarial semanal em função do grau de instrução",
titlefont=12, legend=false)
N = 8
data_x = collect(1:N)
X = [ones(N) data_x]
β = X\data_y
s_q = √(sum(abs2, data_y .- β[1] .- β[2] * data_x)/(N-2))
Cov_q = s_q^2 * inv(X' * X)
Var_q = [[1; x] ⋅ (Cov_q * [1; x]) for x in data_x]
8-element Vector{Float64}: 6169.9875992063535 4054.563279478461 2644.280399659866 1939.1389597505686 1939.1389597505686 2644.2803996598677 4054.5632794784633 6169.987599206355
fator = [quantile(TDist(N-2),q) for q = (0.84, 0.975)]
2-element Vector{Float64}: 1.0839756791279644 2.446911851144969
plot(data_y, seriestype = :scatter, xlims=(0,N+1), ylims=(0,2000), color=2,
xticks=0:N+1, xaxis = "nível de instrução", yaxis="salário (USD)",
label="salário", title="Média salarial semanal em função do grau de instrução",
titlefont=12, legend=:topleft)
plot!(data_x,β[1] .+ β[2] * data_x, grid=false, yerror=fator[2]*sqrt.(Var_q),
ribbon=fator[1]*sqrt.(Var_q), fillalpha=0.2, label="modelo", color=3, legend=:topleft)
Obtenha $\operatorname{Var}(\boldsymbol\beta)$ explicitamente em termos de $N$ e de $(x_i)_{i=1}^N$.
No caso de apenas duas amostras, com $x_2 - x_1 = d$ denotando a distância entre os dois pontos/momentos/condições de medição, escreva $\operatorname{Var}(\boldsymbol\beta)$ explicitamente em função de $d$ e $x_1$, observe a influência da distância $d$ na variância e encontre os limites dessa variância quanto $d\searrow 0$ e $d\nearrow \infty$.
Refaça os exercícios do caderno "Modelos redutíveis ao caso linear nos parâmetros e aplicações", calculando e exibindo os intervalos de confiança em relação às variáveis usadas para linearizar os problemas.
Morris H. DeGroot, Mark J Schervish, "Probability and Statistics", Pearson Education 2012.
John R. Taylor, An Introduction to Error Analysis. The Study of Uncertainties in Physical Measurements. University Science Books, 1997.