0001 function varargout = fpolyfit(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 if nargin == 2 && ischar(varargin{2}) && strcmp(varargin{2}, 'INFO')
0025 varargout{1} = getInfo();
0026 return
0027 end
0028
0029
0030 in_names = cell(size(varargin));
0031 for ii = 1:nargin,in_names{ii} = inputname(ii);end
0032
0033
0034 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0035 [pl, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0036
0037
0038 pl = combine(pl, getDefaultPlist);
0039
0040
0041 orders = find(pl, 'powers');
0042
0043
0044 for j=1:numel(as)
0045 if ~isa(as(j).data, 'fsdata')
0046 warning('!!! Can only fit to fsdata objects. Skipping AO %s', ao_invars{j});
0047 as(j) = [];
0048 else
0049
0050 x = as(j).data.getX;
0051
0052 if ~isempty(N)
0053 [p,s] = polyfit(x, as(j).data.y,N);
0054 else
0055 p = coeffs;
0056 s = [];
0057 end
0058
0059 y = polyval(p, as(j).data.getX);
0060
0061 as(j).data = as(j).data.setXY(x, y);
0062
0063 as(j) = as(j).setName(sprintf('polyfit(%s)', ao_invars{j}));
0064
0065 as(j).procinfo = plist('coeffs', p, 'S', s);
0066
0067 as(j) = as(j).addHistory(getInfo, pl, ao_invars(j), as(j).hist);
0068 end
0069 end
0070
0071
0072 varargout{1} = as;
0073
0074
0075
0076
0077
0078 function [fit_out,dfit_out,C,chi2,N_DOF]=polynom_fit(x,y,orders,varargin)
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114 VERSION = '$Id: polynom_fit.m,v 1.3 2008/03/20 10:25:28 dan Exp $';
0115
0116 plt=1;
0117 talk=1;
0118 err_given=0;
0119
0120 xerr=0;
0121 center=0;
0122 first_time=1;
0123
0124 last_time=1;
0125 sigtest = 0;
0126 logplot = 0;
0127
0128 if ~isempty(varargin)
0129 j = 1;
0130 while j<=length(varargin)
0131 switch lower(repr(varargin{j}))
0132 case 'xran'
0133 xran = varargin{j+1};
0134 j = j+1;
0135 case 'err'
0136 err = varargin{j+1};
0137 err_given = 1;
0138 j = j+1;
0139 case 'errx'
0140 xerr = 1;
0141 errx = varargin{j+1};
0142 last_time=0;
0143 j = j+1;
0144 case 'nopl'
0145 plt=0;
0146 case 'nobs'
0147 talk=0;
0148 case 'center'
0149 center=1;
0150 case 'one_more_time'
0151 first_time=0;
0152 case 'sigtest'
0153 sigtest = 1;
0154 case 'logplot'
0155 logplot = 1;
0156 otherwise
0157 warning('Unknown option "%s"',repr(varargin{j}));
0158 end
0159 j=j+1;
0160 end
0161 end
0162
0163 if exist('xran')
0164 in= x>=xran(1) & x<=xran(2);
0165 xin=x(in);
0166 yin=y(in);
0167 if talk
0168 disp([' Fitting data between x = ' num2str(xran(1)) ' and ' num2str(xran(2)) ')']);
0169 end
0170 else
0171 in=ones(size(x));
0172 xin=x;
0173 yin=y;
0174 if talk
0175 disp([' Fitting all data (from x = ' num2str(min(x)) ' to ' num2str(max(x)) ')']);
0176 end
0177 end
0178
0179 if ~exist('err')
0180 err=ones(length(xin),1);
0181 if talk
0182 disp(' Assuming all points have unity error');
0183 end
0184 elseif length(err)==1
0185 if talk
0186 disp([' Assuming all points have error dy = ' num2str(err)]);
0187 end
0188 err=err*ones(length(xin),1);
0189 else
0190 if talk
0191 disp(' Errors entered point by point');
0192 end
0193 if length(xin)<length(err)
0194 err=err(in);
0195 end
0196 end
0197
0198 if center
0199 x_c=mean(xin);
0200 else
0201 x_c=0;
0202 end
0203
0204
0205 F=[];
0206 for j=1:length(orders)
0207 f=(xin-x_c).^orders(j);
0208 F=[F f];
0209 end
0210
0211 M=size(F,2);
0212 G=[];
0213 V=[];
0214 for i=1:M
0215 V(i) = sum (F(:,i) .* yin ./ err.^2);
0216 for j=1:M
0217 G(i,j)=sum(F(:,i) .* F(:,j) ./ err.^2);
0218 end
0219 end
0220 V=V';
0221
0222
0223
0224
0225 C=inv(G);
0226 fit_out= C * V;
0227
0228 y_fit= F*fit_out;
0229 dy_res=yin-y_fit;
0230 dy_mean=sum(dy_res)/(length(yin)-length(fit_out));
0231
0232
0233 N_DOF=length(yin) - length(fit_out);
0234 chi2=sum(dy_res.^2 ./ err.^2) / N_DOF;
0235
0236 if (err_given==0)
0237
0238
0239 C=C*chi2;
0240 end
0241
0242 dfit_out=[];
0243 for i=1:length(fit_out)
0244 dfit_out(i)=sqrt(C(i,i));
0245 end
0246 dfit_out=dfit_out';
0247
0248
0249 sigma=sqrt(sum(dy_res.^2)/N_DOF);
0250
0251 if talk && first_time
0252 disp(' Fit model: ');
0253 if x_c==0
0254 polynomarg='x';
0255 else
0256 polynomarg='(x - x_c)';
0257 end
0258
0259 for j=1:length(orders)
0260 if orders(j)==0
0261 coeff='a0';
0262 term=coeff;
0263 elseif orders(j)<0
0264 expon=['(' int2str(orders(j)) ')'];
0265 coeff=['am' int2str(abs(orders(j))) ];
0266 term=[coeff ' * ' polynomarg '^' expon];
0267 else
0268 expon=int2str(orders(j));
0269 coeff=['a' int2str(orders(j))];
0270 term=[coeff ' * ' polynomarg '^' expon];
0271 end
0272 if j>1
0273 model=[model ' + ' term];
0274 else
0275 model=term;
0276 end
0277 end
0278 disp([' ' model]);
0279 if center
0280 disp([' using x_c = ' num2str(x_c)]);
0281 end
0282 end
0283 if talk
0284 disp(' ');
0285 disp(' Fit results: ');
0286
0287 if err_given
0288 disp([' Fit chi^2 = ' num2str(chi2) ' per DOF (' int2str(N_DOF) ' DOF)']);
0289 else
0290 disp(' Assume good fit for fit parameter uncertainties');
0291 end
0292 disp([' Average fit residual = ' num2str(dy_mean) ]);
0293 disp([' RMS residual = ' num2str(sigma)]);
0294 disp(' Fit parameters: ');
0295 if length(orders)>=0
0296 for j=1:length(orders)
0297 disp([' Polynomial order ' int2str(orders(j)) ': ' num2str(fit_out(j)) ' +/- ' ...
0298 num2str(dfit_out(j))]);
0299 end
0300 end
0301 disp(' ');
0302 end
0303
0304 if plt && last_time
0305 figure
0306 axtop=axes('position',[.15 .5 .75 .4],'fontsize',14);
0307 if err_given
0308 errorbar(xin,yin,err,'b.');
0309 else
0310 plot(xin,yin,'b.');
0311 end
0312 if logplot
0313 set(gca,'XScale','log')
0314 set(gca,'YScale','log')
0315 end
0316 N=1000;
0317 x_mod=min(xin) + (max(xin)-min(xin))/N*(0:N)';
0318 y_mod=model_eval(orders,fit_out,x_mod);
0319 hold on
0320
0321 plot (x_mod,y_mod,'r');
0322 grid on;
0323 axbot=axes('position',[.15 .1 .75 .3],'fontsize',14);
0324 plot (xin,dy_res,'r');
0325 if err_given
0326 errorbar(xin,dy_res,err,'b.');
0327 else
0328 plot(xin,dy_res,'b.');
0329 end
0330 if logplot
0331 set(gca,'XScale','log')
0332 set(gca,'YScale','log')
0333 end
0334 grid on;
0335 end
0336
0337 if talk
0338 disp(' ');
0339 end
0340
0341 if xerr
0342
0343 disp(' ');
0344 if talk
0345 disp(' Second run to account for errors in the x variable');
0346 if length(errx)==1
0347 disp([' Applying error dx = ' num2str(errx) ' to all points']);
0348 end
0349 end
0350 if (err_given==0)
0351 disp(' Strange choice, considering x errors without the y errors...');
0352 end
0353 dy_dx_mod=zeros(length(x),1);
0354 for j=1:length(orders)
0355 dy_dx_mod=dy_dx_mod + orders(j)* fit_out(j) * xin.^(orders(j)-1);
0356 end
0357 if length(errx)==1
0358 errx=ones(size(xin))*errx;
0359 end
0360 erry_x = sqrt(dy_dx_mod.^2 .* errx.^2);
0361 errs=sqrt(err.^2 + erry_x.^2);
0362
0363 args = {'one_more_time'};
0364 if ~plt
0365 args = [args 'nopl'];
0366 end
0367 if ~talk
0368 args = [args 'nobs'];
0369 end
0370 if center
0371 args = [args 'center'];
0372 end
0373 [fit_out,dfit_out,C,chi2,N_DOF]=polynom_fit(xin,yin,orders,'err',errs,args{:});
0374 end
0375
0376 if sigtest && length(orders)>1
0377 sumsqt = sum((yin-y_fit).^2./err.^2);
0378 sigIdx = zeros(length(orders),1);
0379
0380 if plt
0381 fig = figure;
0382 plot(x,y,'*','Color',0.5*[1,1,1])
0383 hold all
0384 plot(xin,yin,'*','Color',0*[1,1,1])
0385 if logplot
0386 set(gca,'XScale','log')
0387 set(gca,'YScale','log')
0388 end
0389 grid on
0390 plot(xin,y_fit)
0391 legentries{1} = 'Original data (all)';
0392 legentries{2} = 'Original data (selected for fit)';
0393 legentries{3} = ['Complete fit (orders: ' num2str(orders) ')'];
0394 end
0395
0396 args = {};
0397 for j=1:length(varargin)
0398 if ~strcmpi(varargin{j},'sigtest')
0399 args = [args varargin{j}];
0400 end
0401 end
0402
0403 for j=1:length(orders)
0404 disp(' ')
0405 disp(['sum = ' num2str(sumsqt)])
0406 disp(' ')
0407 disp(['Analizing order ' num2str(orders(j))]);
0408 ordersi = orders(orders~=orders(j));
0409 [fiti,dfiti,Ci,chii,N_DOFi] = polynom_fit(x,y,ordersi,'nopl','nobs',args{:});
0410 yfiti = model_eval(ordersi,fiti,xin);
0411 sumsq(j) = sum((yin-yfiti).^2./err.^2);
0412 f(j) = ((sumsq(j)-sumsqt)/(N_DOFi-N_DOF))/(sumsqt/N_DOF);
0413 p(j) = 1 - fcdf(f(j),(N_DOFi-N_DOF),N_DOF);
0414
0415 disp(['sum = ' num2str(sumsq(j))])
0416 disp(['f = ' num2str(f(j))])
0417 disp(['p = ' num2str(p(j))])
0418 if p(j) < 0.05
0419 sigIdx(j) = 1;
0420 disp('SIGNIFICATIVE')
0421 else
0422 disp('NOT SIGNIFICATIVE')
0423 end
0424 if plt
0425 figure(fig)
0426 plot(xin,yfiti)
0427 legentries{j+3} = ['Fit without order ' num2str(orders(j))];
0428 end
0429 end
0430 sigIdx = logical(sigIdx);
0431
0432 [fite,dfite,Ce,chie,N_DOFe] = polynom_fit(x,y,orders(sigIdx),'nopl','nobs',args{:});
0433 if plt
0434 yfite = model_eval(orders(sigIdx),fite,xin);
0435 figure(fig)
0436 plot(xin,yfite,'linewidth',2)
0437 legentries = [legentries ['Fit with relevant orders ( ' num2str(orders(sigIdx)) ')']];
0438 legend(legentries)
0439 end
0440
0441 disp('Summary:')
0442 k = 1;
0443 for j=1:length(orders)
0444 if sigIdx(j)
0445 fit_out(j) = fite(k);
0446 dfit_out(j) = fite(k);
0447 k = k+1;
0448 disp(['Order ' num2str(orders(j)) ': SIGNIFICATIVE'])
0449 else
0450 fit_out(j) = NaN;
0451 dfit_out(j) = NaN;
0452 disp(['Order ' num2str(orders(j)) ': NOT SIGNIFICATIVE'])
0453 end
0454 end
0455 C = Ce;
0456 chi2 = chie;
0457 N_DOF = N_DOFe;
0458 end
0459 end
0460
0461 end
0462
0463