Home > classes > @ao > fpolyfit.m

fpolyfit

PURPOSE ^

FPOLYFIT polynomial fit of powers of f.

SYNOPSIS ^

function varargout = fpolyfit(varargin)

DESCRIPTION ^

 FPOLYFIT polynomial fit of powers of f.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: FPOLYFIT polynomial fit of powers of f.

 CALL:        b = fpolyfit(a, pl)

 PARAMETERS:  'powers' - vector of powers of f to fit

 The following call returns an minfo object about this method:

 >> info = ao.getInfo('fpolyfit')

 VERSION:     $Id: polyfit.m,v 1.12 2008/06/13 10:48:14 hewitson Exp $
 
 HISTORY: 19-06-08 M Hewitson
             Creation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = fpolyfit(varargin)
0002 % FPOLYFIT polynomial fit of powers of f.
0003 %
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %
0006 % DESCRIPTION: FPOLYFIT polynomial fit of powers of f.
0007 %
0008 % CALL:        b = fpolyfit(a, pl)
0009 %
0010 % PARAMETERS:  'powers' - vector of powers of f to fit
0011 %
0012 % The following call returns an minfo object about this method:
0013 %
0014 % >> info = ao.getInfo('fpolyfit')
0015 %
0016 % VERSION:     $Id: polyfit.m,v 1.12 2008/06/13 10:48:14 hewitson Exp $
0017 %
0018 % HISTORY: 19-06-08 M Hewitson
0019 %             Creation
0020 %
0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0022 
0023 %% Check if this is a call for parameters
0024 if nargin == 2 && ischar(varargin{2}) && strcmp(varargin{2}, 'INFO')
0025   varargout{1} = getInfo();
0026   return
0027 end
0028 
0029 % Collect input variable names
0030 in_names = cell(size(varargin));
0031 for ii = 1:nargin,in_names{ii} = inputname(ii);end
0032 
0033 % Collect all AOs and plists
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 % Combine plists
0038 pl = combine(pl, getDefaultPlist);
0039 
0040 % get orders
0041 orders = find(pl, 'powers');
0042 
0043 % Loop over AOs
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     % Get x values
0050     x = as(j).data.getX;
0051     % Fit polynomial
0052     if ~isempty(N)
0053       [p,s] = polyfit(x, as(j).data.y,N);
0054     else
0055       p = coeffs;
0056       s = [];
0057     end
0058     % make new values from these coefficients
0059     y = polyval(p, as(j).data.getX);
0060     % Set data
0061     as(j).data = as(j).data.setXY(x, y);
0062     % Set output name
0063     as(j) = as(j).setName(sprintf('polyfit(%s)', ao_invars{j}));
0064     % Set extra processing info
0065     as(j).procinfo = plist('coeffs', p, 'S', s);
0066     % Add history
0067     as(j) = as(j).addHistory(getInfo, pl, ao_invars(j), as(j).hist);
0068   end
0069 end
0070 
0071 % Set outputs
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     % [fit_out,dfit_out,C,chi2,N_DOF]=polynom_fit(x,y,orders,varargin)
0080     %
0081     % LSQ routine that fits data x(t) to a polynomial with terms of order
0082     % given by vector "orders"
0083     %
0084     % y(x) = a_ord1*x^ord_1 + a_ord2*x^ord_2 +  ...
0085     %       for instance [0,1,-2] fits to a0 + a1*x + a2./x.^2
0086     %
0087     % output:   fitout (fit parameters): [a_ord1; a_ord2; ... ...]
0088     %           NB: the order of the coefficients corresponds to the order in
0089     %           which they were listed in the "orders" vector (confusing!)
0090     %
0091     %           dfit_out (fit uncertainties):   [da_ord1; da_ord2; ... ]
0092     %           C   full covariance matrix
0093     %           chi2    fit chi^2 (meaningless if fit errors not given)
0094     %           N_DOF number of degrees of freedom
0095     %
0096     % options:
0097     % 'err'     specify data uncertainties (next arg either a vector of errors
0098     %           or a single number to be applied to all points ... otherwise
0099     %           unity errors and a good fit are assumed)
0100     % 'errx'    specify data uncertainties in x variable (next argument
0101     %           contains x errors (column vector or a single number)
0102     % 'xran'    next arg is data range (fit only data in the time range)
0103     % 'center'  fits polynomials to the function a_ordj*(x-x_center)^ordj, where x_center
0104     %           is the mean of the in-range x values
0105     % 'sigtest' performs an F-test to evaluate fit coefficients significance with
0106     %           5% threshold, and set to 0 non-significative coefficients.
0107     %           ATTENTION: this feature is still experimental, and may not make
0108     %           sense because we are not using orthogonal polynoms to fit data
0109     % 'logplot' asks for logaritmic plots
0110     % 'nopl', 'nobs' turn off the plotting and text summary outputs
0111     %
0112     % bw 14/12/2006
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; % assume errors are NOT specified, and
0119     % thus we assign unity error to each point
0120     xerr=0; % assume that there are no x errors
0121     center=0;
0122     first_time=1; % internal variable to turn off certain things the second time round
0123     % when evaluating the x errors
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     % construct matrix of functions:
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=[]; % data vector
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'; % it automatically makes V a row vector, we want a column vector
0221 
0222     % "function" matrix
0223     % "data" vector:
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     % assign fit uncertainties, if necessary assuming a good fit (chi^2 = 1)
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       % here we assume a good fit, and scale the covariance matrix such that chi^2 would be = 1
0238       % this allows us to get an idea of the fit errors
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';  % make it a column vector
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       % use model (obtained from first fit) to convert x errors into y errors
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 % End polynom_fit
0460 
0461 end % End fpolyfit
0462 
0463 % END

Generated on Mon 25-Aug-2008 19:02:29 by m2html © 2003