mc

a = 2;
b = 3;
x = -1:0.01:1;
xa = ao(plist('dtype', 'cdata', 'yvals', x));
  
n = 1*randn(size(x));
sn = 2;
na = sn.*ao(plist('dtype', 'cdata', 'yvals', n));

y0 = a.*(xa + b).^3;
ya = y0 + na;

%% linearisation


% ylin = dyda.*da + dydb.*db + y0;

da = 0.2;
db = -0.3;
ag = a+da;
bg = b+db;
y0 = ag.*(xa + bg).^3;

dyda = (xa+bg).^3;
dydb = ag.*(xa+bg).^2;

tgtData = ya-y0;
fitData = lscov(dyda, dydb, tgtData);

dfit = fitData - tgtData;
iplot(dfit)
iplot(hist(dfit))

params = find(fitData.procinfo, 'P')
af = ag + params(1)
bf = bg + params(2)

errs = find(fitData.procinfo, 'STD');

afe = 100*errs(1)./af
bfe = 100*errs(2)./bf

%% Non-linear fit

xy = ao(plist('xvals', x, 'yvals', ya.y));
pl = plist('Function', 'P(1).*(Xdata + P(2)).^3', 'P0', [ag bg]);
nlfit = curvefit(xy, pl);