
def resid(pars, data):
    model = pars.slope.value * data.x + pars.intercept.value
    pars.iter = pars.iter+1
    return  model - data.y
#enddef

print( 'Fits to line')
mdat = group()
mdat.x = arange(100)
mdat.y = 21.3556 - 0.122*mdat.x + 6.e-4*(mdat.x-40)**2
mdat.y = mdat.y + random.normal(size=len(mdat.x), scale=1.75)


######
print( '# Version 1: use minimize() -- most general, overkill for linear fit')
params = param_group(slope = guess(-1),  intercept = guess(0))
params.iter = 0
model = mdat.y + resid(params, mdat)
newplot(mdat.x, mdat.y)
plot(mdat.x, model)

#
out = minimize(resid, params, args=(mdat,))

# newplot(mdat.x, mdat.y)
model = resid(params, mdat) + mdat.y
newplot(mdat.x, mdat.y)
plot(mdat.x, model)
print( fit_report(out))

######
print( '# Version 2: use fitpeak() -- simple wrapper for minimize(), still overkill')
myfit = fit_peak(mdat.x, mdat.y, 'linear')

plot(myfit.x, myfit.y, marker='+', label='data',
     xlabel='x', ylabel='y', show_legend=True, new=True)
plot(myfit.x, myfit.fit_init, label='init')
plot(myfit.x, myfit.fit, label='best fit')

print( fit_report(myfit, min_correl=0.3))

######
print( '# Version 3: use linregress() for linear regression!')
slope, intercept, r_val, p_val, err = linregress(mdat.x, mdat.y)
print( 'Slope=%.6f, Intercept=%.6f' % ( slope, intercept))
plot(mdat.x, intercept + slope*mdat.x, label='regression',
     style='dashed', color='red')
