3.3.4 Aproximação de Funções — Regressões — Regressão Não-Linear
1. Introdução
Nos posts anteriores, vimos como ajustar algumas funções simples pelo método de mínimos quadrados. Antes de apresentarmos o método para regressão não-linear, estudaremos a formulação geral do método de mínimos quadrados, o que nos permite ajustar pontos para funções lineares e não-lineares.
1.1. Mínimos Quadrados Linear Geral
Todos os métodos de regressão apresentados nos posts anteriores — linear simples, linear múltiplo e polinomial — pertencem à um mesmo modelo geral de mínimos quadrados, descrito por
onde cada é uma função-base. Como ilustração, temos o caso em que , , , , é exatamente a função usada no método de regressão linear. Além disso, temos o caso em que cada função-base é um monômio, tal como , , , , , , que constitui o caso de regressão polinomial.
Matricialmente, a equação 1 pode ser expressa como
$$ Y = ZA + E $$
onde é a matriz dos valores calculados das funções-base nos valores medidos de cada variável independente. Isto é, cada variável corresponde à uma linha dessa matriz e cada coluna é uma amostragem. Assim,
$$
Z =
\left[ \begin{array}{cccc}
z_{01} & z_{11} & \cdots & z_{m1} \\
z_{02} & z_{12} & \cdots & z_{m2} \\
\vdots & \vdots & \ddots & \vdots \\
z_{0n} & z_{1n} & \cdots & z_{mn}
\end{array} \right]
$$
onde é o número de variáveis do modelo e é o número de pontos amostrados. O vetor contém os valores os valores observados das variáveis independentes
$$
Y =
\left[ \begin{array}{c}
y_{1} \\
y_{2} \\
\vdots \\
y_{n}
\end{array} \right]
$$
O vetor contém os coeficientes que desejamos ajustar no modelo
$$
A =
\left[ \begin{array}{c}
a_{0} \\
a_{1} \\
\vdots \\
a_{m}
\end{array} \right]
$$
e é o vetor de resíduos (erros)
$$
E =
\left[ \begin{array}{c}
\epsilon_{1} \\
\epsilon_{2} \\
\vdots \\
\epsilon_{n}
\end{array} \right]
$$
Assim, como demonstrado nos posts anteriores, o objetivo da regressão é minimizar o erro, que fazemos através da minimização da seguinte soma dos quadrados dos resíduos:
$$ S_r = \sum_{i=1}^{n} \left(y_i – \sum_{j=0}^{m} a_j z_{ji} \right)^2 $$
Essa soma pode ser minimizada tomando-se suas derivadas parciais em relação a cada um dos coeficientes a serem ajustados e igualando-se a equação resultante a zero. Desta maneira, temos que os coeficientes são obtidos ao resolver o sistema
$$ Z^T Z A = Z^T Y $$
2. Regressão Não-linear
Como nos casos dos mínimos quadrados linear, a regressão não-linear baseia-se na obtenção dos coeficientes que minimizam os resíduos, com a diferença que o processo de minimização é feito iterativamente por meio do método de Gauss-Newton. Esse método consiste em expressar a função não-linear através de uma aproximação linear por expansão em série de Taylor. Assim, a teoria dos mínimos quadrados pode ser usada para obter-se novas estimativas dos parâmetros desconhecidos.
Para ilustrar esse método, a relação entre a equação não-linear e os dados pode ser expressa como
$$ y_i = f(x_i; a_0, a_1, a_2, \ldots, a_m) + \epsilon_i $$
onde é uma amostra da variável dependente; é a função independente da variável e dependente dos parâmetros ; e é um resíduo qualquer. Por conveniência, omitimos os parâmetros,
$$ y_i = f(x_i) + \epsilon_i $$
Expandimos esse modelo por uma série de Taylor em torno dos parâmetros e truncamos após a primeira derivada. Isto é, para um modelo com parâmetros, temos
$$
f(x_i)_{j+1} = f(x_i)_j +
\frac{\partial f(x_i)_j}{\partial a_0} \Delta a_0 +
\frac{\partial f(x_i)_j}{\partial a_1} \Delta a_1 +
\cdots +
\frac{\partial f(x_i)_j}{\partial a_m} \Delta a_m
$$
onde é a aproximação inicial, é a previsão, e . Com isso, linearizamos o
modelo inicial com relação aos parâmetros. Levando a Equação 4 na Equação 5, temos
$$
y_i – f(x_i)_j =
\frac{\partial f(x_i)_j}{\partial a_0} \Delta a_0 +
\frac{\partial f(x_i)_j}{\partial a_1} \Delta a_1 +
\cdots +
\frac{\partial f(x_i)_j}{\partial a_m} \Delta a_m +
\epsilon_i
$$
ou, na forma matricial, temos um sistema igual da Equação 2,
$$ D = Z A + E $$
onde é a matriz das derivadas parciais da função calculada na aproximação inicial,
$$
Z =
\left[ \begin{array}{cccc}
\frac{\partial f_1}{\partial a_0} &
\frac{\partial f_1}{\partial a_1} &
\cdots &
\frac{\partial f_1}{\partial a_m} \\
\frac{\partial f_2}{\partial a_0} &
\frac{\partial f_2}{\partial a_1} &
\cdots &
\frac{\partial f_2}{\partial a_m} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial f_n}{\partial a_0} &
\frac{\partial f_n}{\partial a_1} &
\cdots &
\frac{\partial f_n}{\partial a_m} \\
\end{array} \right]
$$
onde é o número de amostrars e é a derivada parcial da função com relação ao k-ésimo parâmetro calculada no i-ésimo ponto dado. O vetor contém as diferenças entre as medidas e os valores da função ajustada,
$$
D =
\left[ \begin{array}{c}
y_1 – f(x_1) \\
y_2 – f(x_2) \\
\vdots \\
y_n – f(x_n)
\end{array} \right]
$$
e o vetor contém a variação nos valores dos parâmetros que estão sendo ajustados,
$$
A =
\left[ \begin{array}{c}
\Delta a_0 \\
\Delta a_1 \\
\vdots \\
\Delta a_n \\
\end{array} \right]
$$
Assim, a Equação 7 pode ser reescrita como a seguinte equação normal (conforme apresentado na Equação normal:
$$ Z^T Z A = Z^T D $$
A resolução deste sistema ajusta o conjunto de coeficientes da função não-linear.
Devemos notar que o vetor dos coeficientes ( ) resolve uma variação dos coeficientes e não os valores dos coeficientes em si. Isso descreve um processo iterativo, ou seja, esse conjunto de equações pode ser resolvidas sucessivas vezes, substituindo os valores de pelos novos valores melhorados como em
| $$ a_{0,j+1} = a_{0,j} + \Delta a_0 $$ |
| $$ a_{0,j+1} = a_{1,j} + \Delta a_1 $$ |
| $$ \vdots $$ |
| $$ a_{m,j+1} = a_{m,j} + \Delta a_m $$ |
</center>
Este processo iterativo é repetido até que a solução convirja à um valor de erro tolerado, que pode ser calculado utilizando a seguinte equação
$$ \left|\epsilon_a\right|_k = \left|\frac{a_{k,j+1} – a_{k,j}{a_{k,j+1}\right| $$
3. Implementação
Para ilustrarmos a utilização da regressão linear, implementaremos a seguinte função
$$ f(x) = \alpha_0(1 – e^{-\alpha_1 x}) $$
ilustrando cada passo e o respectivo código, seguindo a ordem do algoritmo descrito na seção anterior. O primeiro passo é definir uma aproximação inicial para $$\alpha_0 $$ e $$\alpha_1 $$ e selecionar um erro de tolerância $$\epsilon $$ . Como a função possui dos coeficientes a serem ajustados, o vetor dos coeficientes terá apenas dois elementos, ou seja $$A = [\alpha_0, \alpha_1]^T $$ . Além disso, escolhemos $$\epsilon=0.1 $$ .
O próximo passo é obter as derivadas parciais da função com relação aos parâmetros $$\alpha_0 $$ e $$\alpha_1 $$ . Isto é,
$$ \frac{\partial f}{\partial \alpha_0} = 1 – e^{-\alpha_1 x} $$ e $$ \frac{\partial f}{\partial \alpha_1} = \alpha_0 x e^{-\alpha_1 x} $$
Na implementação temos que as funções necessárias $$f $$ , $$\frac{\partial f}{\partial \alpha_0}$$ e $$\frac{\partial f}{\partial \alpha_1}$$ são, respectivamente,
def func(x, a0, a1):
return a0 * (1 - exp(-a1 * x))
def diff_a0(x, a0, a1):
return 1.0 - exp(-a1 * x)
def diff_a1(x, a0, a1):
return a0 * x * exp(-a1 * x)
</div>
As equações acima são usadas para calcular a matriz $$Z $$ , conforme mostrado na seção anterior
$$
Z =
\left[ \begin{array}{cc}
\frac{\partial f_1}{\partial a_0} &
\frac{\partial f_1}{\partial a_1} \\
\frac{\partial f_2}{\partial a_0} &
\frac{\partial f_2}{\partial a_1} \\
\vdots & \vdots \\
\frac{\partial f_n}{\partial a_0} &
\frac{\partial f_n}{\partial a_1} \\
\end{array} \right]
$$
Essa matriz é construída com a seguinte implementação
def get_Z(xpoints, a0=1.0, a1=1.0):
lines = len(xpoints)
Z = [[0 for i in xrange(2)] for j in xrange(lines)]
for i in xrange(lines):
Z[i][0] = diff_a0(xpoints[i], a0, a1)
Z[i][1] = diff_a1(xpoints[i], a0, a1)
return Z
</div>
O passo seguinte consiste em obter o vetor $$D $$ que, conforme demonstrado, será a diferença entre a função ajustada e o valor amostrado
def get_D(y, x, a0, a1):
n = len(y)
D = [0.0 for i in xrange(n)]
for i in xrange(n):
D[i] = y[i] - func(x[i], a0, a1)
return D
</div>
Agora, de posse de $$Z $$ e $$D $$ , podemos obter o vetor $$A $$ resolvendo o sistema da Equação sys. Os passos para resolução desse sistema são listados nos comandos abaixo
def calcule_alphas(y, x, alphas):
Z = get_Z(x, alphas[0], alphas[1])
Zt = transpose(Z)
ZtZ = matmul(Zt, Z)
D = get_D(y, x, alphas[0], alphas[1])
ZtD = matmul(Zt, D)
A = cholesky(ZtZ, ZtD)
return [i + j for (i, j) in zip(alphas, A)]
</div>
Logo, os coeficientes já estão mais ajustados em relação à aproximação inicial. Contudo, esse procedimento pode ser repetido para melhorar ainda mais os parâmetros $$\alpha $$ . Logo, definimos uma nova função que ajusta tais parâmetros através de sucessivas substituições de $$\alpha $$:
def iterative_general_least_squares(x, y):
alphas = [1.0, 1.0]
iterCount = 0
errno = 1
while errno > errto and iterCount < imax:
memalpha = alphas
alphas = calcule_alphas(y, x, alphas)
iterCount += 1
errno = calcule_error(alphas, memalpha)
return alphas
</div>
Onde a variável memalpha armazena o $$\alpha $$ da iteração anterior para o calculo do erro; itercount é um contador usado para limitar o número máximo de repetições, caso haja uma convergência muito lenta do algoritmo.
Um exemplo completo que utiliza essas funções pode ser encontrado em http://www.sawp.com.br/code/regression/nonlinear_regression.py.
4. Copyright
Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.