# manage data and fit
import pandas as pd
import numpy as np
# scipy ODR : orthogonal distance regression
from scipy.odr import ODR, Model, Data, RealData
# style and notebook integration of the plots
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set(style="whitegrid")
def func(beta, x):
y = beta[0]+beta[1]*x+beta[2]*x**3
return y
npts = 100
np.random.seed(160810)
x = np.linspace(-3, 2, npts)
y = func([-2.3, 7.0, -4.0], x)
sigma_x = .3
sigma_y = 2
# add some noise
x += np.random.normal(scale=sigma_x, size=100)
y += np.random.normal(scale=sigma_y, size=100)
data = RealData(x, y, 0.3, 0.1)
model = Model(func)
odr = ODR(data, model, [1,0,0])
odr.set_job(fit_type=0)
lsq_output = odr.run()
print(" stop reason:", lsq_output.stopreason)
print(" params:", lsq_output.beta)
print(" info:", lsq_output.info)
print("Sum of squares error:", lsq_output.sum_square)
print(" sd_beta:", lsq_output.sd_beta)
print(" sqrt(diag(cov):", np.sqrt(np.diag(lsq_output.cov_beta)))
# check convergence and run again the algorithm if it is not reached
if lsq_output.info != 1:
print("\nRestart ODR till convergence is reached")
i = 1
while lsq_output.info != 1 and i < 100:
print("restart %3d %12.7f" % (i, lsq_output.sum_square))
lsq_output = odr.restart()
i += 1
print(" stop reason:", lsq_output.stopreason)
print(" params:", lsq_output.beta)
print(" info:", lsq_output.info)
print("Sum of squares error:", lsq_output.sum_square)
print(" sd_beta:", lsq_output.sd_beta)
print(" sqrt(diag(cov):", np.sqrt(np.diag(lsq_output.cov_beta)))
stop reason: ['Iteration limit reached'] params: [-2.09368488 6.0257117 -2.54574759] info: 4 Sum of squares error: 5405.197746282522 sd_beta: [ 0.22559351 0.80177874 0.95293199] sqrt(diag(cov): [ 0.03022085 0.1074075 0.12765622] Restart ODR till convergence is reached restart 1 5405.1977463 restart 2 4325.2070339 restart 3 1045.5356553 restart 4 743.3175690 restart 5 710.8233373 restart 6 682.0088933 restart 7 602.7134451 restart 8 578.8585227 restart 9 558.9136790 restart 10 539.9567887 restart 11 522.0568694 restart 12 505.0221792 restart 13 488.8286431 restart 14 473.3871052 restart 15 458.6550721 restart 16 444.5728008 restart 17 431.1303083 restart 18 418.2785524 restart 19 406.0148268 restart 20 394.2378398 restart 21 383.0217229 restart 22 372.3116187 restart 23 362.0182202 restart 24 352.2750168 restart 25 342.8414632 restart 26 333.9524660 restart 27 325.3301467 restart 28 317.2122579 restart 29 309.3065211 restart 30 301.8936369 restart 31 294.6241475 restart 32 287.8565250 restart 33 281.1543201 restart 34 274.9756621 restart 35 250.7690432 restart 36 245.9163778 restart 37 241.1998286 restart 38 236.6336790 restart 39 232.2524033 restart 40 228.0534695 restart 41 224.0110504 restart 42 220.1177572 restart 43 216.3814990 restart 44 212.8300817 restart 45 209.4103793 restart 46 206.1065870 restart 47 202.9507254 restart 48 199.9235477 restart 49 197.0054147 restart 50 194.2163804 restart 51 191.5362613 restart 52 188.9598303 restart 53 186.4947133 restart 54 184.1129010 restart 55 181.8341270 restart 56 179.6356961 restart 57 177.5532439 restart 58 175.5140346 restart 59 173.5668132 restart 60 171.6843892 restart 61 169.8962720 restart 62 168.1662248 restart 63 166.5032540 restart 64 164.8891762 restart 65 163.3493519 restart 66 161.8492549 restart 67 160.3945986 restart 68 157.0580837 restart 69 153.3941930 restart 70 149.0470063 restart 71 148.1832034 restart 72 146.9758074 restart 73 145.7639864 restart 74 145.2494874 restart 75 145.1131538 restart 76 145.0754139 restart 77 145.0461410 restart 78 144.9885690 restart 79 144.9623475 restart 80 144.9537335 restart 81 144.9443404 restart 82 144.9438557 restart 83 144.9434179 restart 84 144.9430670 restart 85 144.9430219 restart 86 144.9429878 restart 87 144.9429137 stop reason: ['Sum of squares convergence'] params: [ -2.01039971 12.55629008 -6.06416408] info: 1 Sum of squares error: 144.94291057820132 sd_beta: [ 0.48021907 0.89693549 0.49629295] sqrt(diag(cov): [ 0.39285002 0.7337508 0.40599949]
lsq_output.pprint()
Beta: [ -2.01039971 12.55629008 -6.06416408] Beta Std Error: [ 0.48021907 0.89693549 0.49629295] Beta Covariance: [[ 0.15433114 0.15842579 -0.04105628] [ 0.15842579 0.53839023 -0.21754393] [-0.04105628 -0.21754393 0.16483558]] Residual Variance: 1.4942568100845495 Inverse Condition #: 0.037046778050718744 Reason(s) for Halting: Sum of squares convergence
odr = ODR(data, model, [1,0,0])
odr.set_job(fit_type=2)
odr_output = odr.run()
print(" stop reason:", odr_output.stopreason)
print(" params:", odr_output.beta)
print(" info:", odr_output.info)
print("Sum of squares error:", odr_output.sum_square)
print(" sd_beta:", odr_output.sd_beta)
print(" sqrt(diag(cov):", np.sqrt(np.diag(odr_output.cov_beta)))
# check convergence and run again the algorithm if it is not reached
if odr_output.info != 1:
print("\nRestart ODR till convergence is reached")
i = 1
while odr_output.info != 1 and i < 100:
print("restart", i)
odr_output = odr.restart()
i += 1
print(" stop reason:", odr_output.stopreason)
print(" params:", odr_output.beta)
print(" info:", odr_output.info)
print("Sum of squares error:", odr_output.sum_square)
print(" sd_beta:", odr_output.sd_beta)
print(" sqrt(diag(cov):", np.sqrt(np.diag(odr_output.cov_beta)))
stop reason: ['Sum of squares convergence'] params: [-1.15111597 4.5357938 -3.34692538] info: 1 Sum of squares error: 821948.608722443 sd_beta: [ 1.00126769 1.25522586 0.23236659] sqrt(diag(cov): [ 0.01087712 0.01363595 0.00252428]
xn = np.linspace(-3, 2, 50)
plt.errorbar(x, y, marker="o", linestyle="", label="data", xerr=sigma_x, yerr=sigma_y, elinewidth=1, capthick=1)
plt.plot(xn, func(lsq_output.beta, xn), label='leastsq')
plt.plot(xn, func(odr_output.beta, xn), label='odr')
plt.legend()
<matplotlib.legend.Legend at 0x117930f28>