import sys
sys.path.insert(0, '/home/leo/src/fatiando-inversion')
from fatiando.inversion.base import Misfit
import matplotlib.pyplot as plt
import numpy as np
Equação da reta: $ y = a x + b$
x = np.linspace(0, 100, 100)
x
y = 2*x + 30
y
%matplotlib inline
plt.plot(x, y, '.k')
from fatiando import utils
y = utils.contaminate(y, 0.05, percent=True)
plt.plot(x, y, '.k')
Determinar a e b que ajuste esses dados.
$y_1 = a x_1 + b$
$y_2 = a x_2 + b$
...
$y_i = a x_i + b$
$\mathbf{y} = \mathbf{A}\mathbf{p}$
$\mathbf{A}^T \mathbf{A} \mathbf{p} = \mathbf{A}^T \mathbf{y}$
class AjusteReta(Misfit):
def __init__(self, x, y):
super(AjusteReta, self).__init__(data=y, nparams=2, islinear=True)
self.x = x
def predicted(self, p):
a, b = p
return a*self.x + b
def jacobian(self, p):
jac = np.empty((self.ndata, self.nparams))
jac[:, 0] = self.x
jac[:, 1] = 1
return jac
solver = AjusteReta(x, y)
solver.fit()
solver.estimate_
solver.predicted()
solver.hessian(None)
solver.gradient(None)
plt.plot(x, y, '.k', x, solver.predicted(), '-r')
solver.config('newton', initial=[0, 0]).fit()
solver.estimate_
solver.config('levmarq', initial=[0, 0]).fit()
solver.estimate_
solver.config('acor', bounds=[0, 10, -1000, 1000]).fit()
solver.estimate_