1 | function f = ezfit(varargin) |
---|
2 | %EZFIT Fit data with arbitrary fitting function |
---|
3 | % EZFIT(FUN) fits the active curve with the function FUN. See below for |
---|
4 | % the syntax of FUN. If FUN is not specified, 'linear' is used. |
---|
5 | % |
---|
6 | % By default, the first curve in the active figure is fitted - see |
---|
7 | % FITPARAM to change this default behavior. To fit another curve, select |
---|
8 | % it before calling EZFIT. If some data are selected by the "Data |
---|
9 | % Brushing" tool (only for Matlab >= 7.6), only those data are fitted. |
---|
10 | % |
---|
11 | % EZFIT(X,Y,FUN) or EZFIT(Y,FUN) fits the data (X,Y) (or Y) using the |
---|
12 | % function FUN. X and Y must be vectors of equal length. If X is not |
---|
13 | % specified, X=[1, 2, 3...] is assumed. |
---|
14 | % |
---|
15 | % EZFIT(X,Y,FUN), where X is a 1-by-N vector and Y is a 2-by-N matrix, |
---|
16 | % also specifies the weights for Y(2,:) (error bars). By default, when |
---|
17 | % Y is a 1-by-N vector, all the weights are 1. |
---|
18 | % |
---|
19 | % Note that EZFIT only computes the coefficients, but does not display the |
---|
20 | % fit. Use SHOWFIT to display the fit. |
---|
21 | % |
---|
22 | % The function string FUN can be: |
---|
23 | % - the name of a predefined fitting function (see below). |
---|
24 | % - the name of a user-defined fitting function (see EDITFIT). |
---|
25 | % - an equation, in the form 'y(x)=...', where 'x' represents the |
---|
26 | % X-data, and all the other variables are parameters to be fitted |
---|
27 | % ('a', 'x_0', 'tau', ...). Example: 'y(x)=a*sin(b*x)'. If the |
---|
28 | % left-hand-side 'y(x)' is not specified, 'x' is taken for the |
---|
29 | % X-Data. All the parameter names are accepted, except Matlab |
---|
30 | % reserved strings ('sin', 'pi', ...) |
---|
31 | % |
---|
32 | % The predefined fitting functions are: |
---|
33 | % - linear y = m * x |
---|
34 | % - affine or poly1 y = a*x + b |
---|
35 | % - poly{n} y = a0 + a1 * x + ... + an * x^n |
---|
36 | % - power y = c*x^n |
---|
37 | % - sin y = a * sin (b * x) |
---|
38 | % - cos y = a * cos (b * x) |
---|
39 | % - exp y = a * exp (b * x) |
---|
40 | % - log y = a * log (b * x) |
---|
41 | % - cngauss y = exp(-x^2/(2*s^2))/(2*pi*s^2)^(1/2) |
---|
42 | % - cfgauss y = a*exp(-x^2/(2*s^2)) |
---|
43 | % - ngauss y = exp(-(x-x0)^2/(2*s^2))/(2*pi*s^2)^(1/2) |
---|
44 | % - gauss y = a*exp(-(x-x0)^2/(2*s^2)) |
---|
45 | % 'ngauss' is a 2-parameters normalized Gaussian, and 'gauss' is a |
---|
46 | % 3-parameters non-normalized (free) Gaussian. 'cngauss' and 'cfgauss' |
---|
47 | % are centered normalized and centered free Gaussian, respectively. |
---|
48 | % |
---|
49 | % EZFIT is based on Matlab's built-in FMINSEARCH function (Nelder-Mead |
---|
50 | % method), which performs an unconstrained nonlinear minimization of |
---|
51 | % the SSR (sum of squared residuals) with respect to the various |
---|
52 | % parameters. |
---|
53 | % |
---|
54 | % The correlation coefficient R is defined as (SSR/(SSE+SSR))^(1/2), where |
---|
55 | % SSR = sum ((y_fit - mean(y)).^2) % sum of squared residuals |
---|
56 | % SSE = sum ((y_fit - y).^2) % sum of squared errors |
---|
57 | % (see http://mathworld.wolfram.com/CorrelationCoefficient.html) |
---|
58 | % |
---|
59 | % Nonlinear minimization requires starting guesses (or starting estimates) |
---|
60 | % for the fit parameters. By default, all the starting guesses are taken |
---|
61 | % as 1, or, when using predefined fits (e.g., exp, gauss, power...), the |
---|
62 | % starting guesses are determined depending on the range of the data to |
---|
63 | % be fitted. However, in most cases, values closer to the expected |
---|
64 | % result should be specified to "help" the convergence. It is sufficient |
---|
65 | % to choose values that have the correct sign and correct order of |
---|
66 | % magnitude, e.g. 0.01, 1, 100... |
---|
67 | % |
---|
68 | % The starting guesses for the parameters may be specified in two ways: |
---|
69 | % - directly in the string FUN, after the fit definition: |
---|
70 | % 'c0 + a*sin(pi*x/lambda); c0=1; a=0.1; lambda=100' |
---|
71 | % ('!' or '$' may also be used instead of ';'). |
---|
72 | % - by specifying them as an additional input argument for EZFIT: |
---|
73 | % EZFIT(x,y,'c0 + a*sin(pi*x/lambda)',[0.1 1 100]); |
---|
74 | % Note that in this case the parameters must be ordered alphabetically. |
---|
75 | % Note that if both methods are used, only the starting guesses given in |
---|
76 | % the string FUN are considered. |
---|
77 | % |
---|
78 | % By default, Y is fitted in linear mode. If you want to fit LOG(Y) |
---|
79 | % instead, you must specify the option 'log' to the string FUN, separated |
---|
80 | % by the symbol ';' or '$' or '!' (eg, FUN='a*x^n;log'). This is |
---|
81 | % specially useful to fit power laws with equally weighted points in a |
---|
82 | % log scale. If nothing specified, the option 'lin' is used. |
---|
83 | % |
---|
84 | % Example: |
---|
85 | % plotsample('power') |
---|
86 | % and compare the results of: |
---|
87 | % ezfit('power;lin') |
---|
88 | % ezfit('power;log') |
---|
89 | % |
---|
90 | % F = EZFIT(...) also returns a structure F having the following fields: |
---|
91 | % - name name of the fit |
---|
92 | % - eq equation of the fit |
---|
93 | % - param cell array of strings: names of the parameters |
---|
94 | % - m values of the coefficients |
---|
95 | % - m0 initial guess for the coefficients |
---|
96 | % - r correlation coefficient R (Pearson's correlation) |
---|
97 | % - fitmode 'lin' (y is fitted) or 'log' (log(y) is fitted) mode |
---|
98 | % |
---|
99 | % This structure F can be further used with SHOWFIT, DISPEQFIT, |
---|
100 | % SHOWEQBOX, MAKEVARFIT and EDITCOEFF. |
---|
101 | % |
---|
102 | % From F, you can get the values of the fitted parameters. If you want to |
---|
103 | % create in the current Matlab workspace the variables associated to |
---|
104 | % these parameters, use MAKEVARFIT (or set the option 'automakevarfit' |
---|
105 | % to 'on' in FITPARAM). |
---|
106 | % |
---|
107 | % Examples: |
---|
108 | % plotsample damposc |
---|
109 | % f = ezfit('u(t) = c + u_a * sin(2*pi*t/T) * exp(-t/tau); T=5; tau=20'); |
---|
110 | % showfit(f); |
---|
111 | % |
---|
112 | % plotsample poly2 |
---|
113 | % [x,y] = pickdata; |
---|
114 | % f = ezfit(x, y, 'z(v) = poly3'); |
---|
115 | % editcoeff(f); |
---|
116 | % |
---|
117 | % plotsample poly2 |
---|
118 | % f = ezfit('beta(z) = poly2'); |
---|
119 | % showfit(f, 'fitcolor', 'red', 'fitlinewidth', 2); |
---|
120 | % |
---|
121 | % See also SHOWFIT, PLOTSAMPLE, DISPEQFIT, EDITCOEFF, |
---|
122 | % FMINSEARCH, MAKEVARFIT, FITPARAM. |
---|
123 | |
---|
124 | |
---|
125 | % F. Moisy, moisy_at_fast.u-psud.fr |
---|
126 | % Revision: 2.51, Date: 2011/10/04 |
---|
127 | % This function is part of the EzyFit Toolbox |
---|
128 | |
---|
129 | |
---|
130 | % History: |
---|
131 | % 2005/05/12: v1.00, first version. |
---|
132 | % 2005/05/20: v1.10, Use 'eval', with generic functions |
---|
133 | % 2005/05/21: v1.11, added the 'lin','log' options. |
---|
134 | % 2005/05/24: v1.12, option 'log' by default for 'power' and 'exp'. |
---|
135 | % 2005/07/27: v1.13, cosmetics. |
---|
136 | % 2005/09/03: v1.14, check arg. |
---|
137 | % 2005/10/07: v1.15, gaussian fits added (ngauss and fgauss, centered/not) |
---|
138 | % 2005/10/20: v1.16, help text changed. |
---|
139 | % 2005/10/31: v1.20, also returns R. cste and poly{n} fits added. Initial |
---|
140 | % guess defined within the fitting function string. The |
---|
141 | % order of the output parameters is changed. |
---|
142 | % 2005/11/05: v1.21, evaluate strings for initial guess in FUN |
---|
143 | % 2005/12/06: v1.22, opens a dialog box if the polynomial order is not |
---|
144 | % specified. |
---|
145 | % 2006/01/13: v1.24, check for negative data in log mode |
---|
146 | % 2006/01/19: v1.25, bug fixed from 1.24 |
---|
147 | % 2006/01/25: v1.26, check the matlab version |
---|
148 | % 2006/02/08: v2.00, new syntax. The output argument is now a fit |
---|
149 | % structure, and the fitting equation string accepts |
---|
150 | % arbitrary parameter names. |
---|
151 | % 2006/02/14: v2.10, lhs 'y(x)=...' now accepted. Now case sensitive |
---|
152 | % 2006/03/07: v2.11, bug fixed from v1.25 |
---|
153 | % 2006/03/09: v2.20, weighted chi-square criterion (ie: error bars accepted |
---|
154 | % for y) - undocumented |
---|
155 | % 2006/09/04: v2.21, '!' and '$' may be used instead of ';' in the FUN |
---|
156 | % string: this allows to pass the argument like: |
---|
157 | % fit power!log |
---|
158 | % 2006/09/28: v2.22, 'x^1' replaced by 'x' for 1st order polynomial |
---|
159 | % 2006/10/18: v2.30, input raw and column vectors accepted. accepts |
---|
160 | % additional input parameters to change the default settings. |
---|
161 | % 2007/04/17: v2.31, help text improved; standard error messages |
---|
162 | % 2007/05/16: v2.40, weigthed fits (v2.20) now documented, with bug fixed. |
---|
163 | % 2007/08/18: v2.41, help text improved. |
---|
164 | % 2007/09/17: v2.50, guess the initial guess m0 for predefined fits |
---|
165 | % 2011/11/04: v2.51, bug fixed for fits in log coordinates |
---|
166 | % 2012/06/28: v2.52, check version number removed (unstable) |
---|
167 | |
---|
168 | |
---|
169 | |
---|
170 | |
---|
171 | % new v1.26, changed v2.11, removed v2.52 |
---|
172 | % if str2double(version('-release'))<14 |
---|
173 | % error('Ezyfit:ezfit:compatibility','EzyFit requires Matlab 7 or higher.'); |
---|
174 | % end |
---|
175 | |
---|
176 | |
---|
177 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
178 | % fit parameters: (new v2.30) |
---|
179 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
180 | |
---|
181 | |
---|
182 | % loads the default fit parameters: |
---|
183 | try |
---|
184 | fp=fitparam; |
---|
185 | catch |
---|
186 | error('Ezyfit:ezfit:fitparamNotFound','''fitparam.m'' file not found.'); |
---|
187 | end |
---|
188 | |
---|
189 | |
---|
190 | % change the default values of the fit parameters according to the |
---|
191 | % additional input arguments: |
---|
192 | for nopt=1:(nargin-1) |
---|
193 | if any(strcmp(varargin{nopt},fieldnames(fp))) % if the option string is one of the fit parameter |
---|
194 | fp.(varargin{nopt}) = varargin{nopt+1}; |
---|
195 | end |
---|
196 | end |
---|
197 | |
---|
198 | |
---|
199 | |
---|
200 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
201 | % input arguments: |
---|
202 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
203 | |
---|
204 | if nargin==0, % if no input argument: do a linear fit of the current curve |
---|
205 | [x,y,h] = pickdata(fp); |
---|
206 | inputstr='linear'; |
---|
207 | elseif nargin==1 |
---|
208 | if ischar(varargin{1}), % EZFIT('a*x+b') |
---|
209 | inputstr=varargin{1}; |
---|
210 | [x,y,h] = pickdata(fp); |
---|
211 | else % EZFIT(y) |
---|
212 | y=varargin{1}; |
---|
213 | x=1:length(y); |
---|
214 | inputstr='linear'; |
---|
215 | end |
---|
216 | else % 2 or more input arguments |
---|
217 | if ~isnumeric(varargin{1}) % EZFIT('fun',...) |
---|
218 | inputstr=varargin{1}; |
---|
219 | [x,y,h] = pickdata(fp); |
---|
220 | if isnumeric(varargin{2}) % EZFIT('fun',m0,...) |
---|
221 | m0=varargin{2}; |
---|
222 | end |
---|
223 | elseif (isnumeric(varargin{1}) && ~isnumeric(varargin{2})) % EZFIT(y,'a*x+b',...) |
---|
224 | [y,inputstr]=deal(varargin{1:2}); |
---|
225 | x=1:length(y); |
---|
226 | if nargin>2 |
---|
227 | if isnumeric(varargin{3}) % EZFIT(y,'a*x+b',m0,...) |
---|
228 | m0=varargin{3}; |
---|
229 | end |
---|
230 | end |
---|
231 | elseif (isnumeric(varargin{1}) && isnumeric(varargin{2})) % EZFIT(x,y,...) |
---|
232 | [x,y]=deal(varargin{1:2}); |
---|
233 | if nargin>2 |
---|
234 | if ischar(varargin{3}) |
---|
235 | inputstr=varargin{3}; % EZFIT(x,y,'fun',...) |
---|
236 | if nargin>3 |
---|
237 | if isnumeric(varargin{4}) % EZFIT(x,y,'fun',m0,...) |
---|
238 | m0=varargin{4}; |
---|
239 | end |
---|
240 | end |
---|
241 | else |
---|
242 | error('Ezyfit:ezfit:syntaxError','Syntax error. 3rd paramater of EZFIT should be a string.'); |
---|
243 | end |
---|
244 | end |
---|
245 | end |
---|
246 | end |
---|
247 | |
---|
248 | |
---|
249 | |
---|
250 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
251 | % some checks about x and y |
---|
252 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
253 | |
---|
254 | |
---|
255 | % turn all input vectors into row (1xN) vectors (new v2.30): |
---|
256 | if size(x,1)>size(x,2), |
---|
257 | x=x'; |
---|
258 | end |
---|
259 | if size(y,1)>size(y,2), |
---|
260 | y=y'; |
---|
261 | end |
---|
262 | |
---|
263 | |
---|
264 | % check for error bars for y (new v2.20, fixed 2.40) |
---|
265 | if size(y,1)>1 |
---|
266 | dy = y(2,:); % second line = error bars (1/weight) |
---|
267 | y = y(1,:); % first line = data |
---|
268 | else |
---|
269 | dy = ones(1,length(y)); % default error bars: 1 |
---|
270 | end |
---|
271 | |
---|
272 | if length(x)~=length(y), |
---|
273 | error('Ezyfit:ezfit:dimagree','X and Y dimensions must agree.'); |
---|
274 | end |
---|
275 | |
---|
276 | |
---|
277 | |
---|
278 | |
---|
279 | |
---|
280 | |
---|
281 | |
---|
282 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
283 | % processing of the string FUN |
---|
284 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
285 | |
---|
286 | % cleans the input string: |
---|
287 | inputstr = strrep(inputstr,' ',''); |
---|
288 | inputstr = strrep(inputstr,'!',';'); |
---|
289 | inputstr = strrep(inputstr,'$',';'); |
---|
290 | inputstr = strrep([inputstr ';'],';;',';'); % ensures that fun is terminated by a ';'. |
---|
291 | |
---|
292 | % the name of the fit is by default the first part of the input string |
---|
293 | p=findstr(inputstr,';'); p=p(1); |
---|
294 | f.name = inputstr(1:(p-1)); |
---|
295 | |
---|
296 | % separates the first part (fitting function itself) and the remainder: |
---|
297 | p = strfind(inputstr,';'); p=p(1); |
---|
298 | fun = inputstr(1:(p-1)); |
---|
299 | remfun = inputstr((p+1):end); |
---|
300 | |
---|
301 | |
---|
302 | % search for predefined fit or user-defined fit |
---|
303 | usepredefinedfit = ''; |
---|
304 | [defaultfit, userfit] = loadfit; |
---|
305 | for i=1:length(defaultfit), |
---|
306 | if strcmp(fun, defaultfit(i).name); |
---|
307 | fun = defaultfit(i).eq; |
---|
308 | usepredefinedfit = defaultfit(i).name; % new v2.50 |
---|
309 | end |
---|
310 | end |
---|
311 | for i=1:length(userfit), |
---|
312 | if strcmp(fun, userfit(i).name), |
---|
313 | fun = userfit(i).eq; |
---|
314 | end |
---|
315 | end |
---|
316 | |
---|
317 | |
---|
318 | % separates again the first part (fitting function itself) and the remainder: |
---|
319 | % (because the predefined/user-defined part may itself contain ';') |
---|
320 | fun = [strrep(fun,' ','') ';' remfun]; |
---|
321 | p=strfind(fun,';'); p=p(1); |
---|
322 | remfun = fun((p+1):end); |
---|
323 | fun = fun(1:(p-1)); |
---|
324 | |
---|
325 | |
---|
326 | % recognize if a lhs is present |
---|
327 | peq = strfind(fun,'='); |
---|
328 | if ~isempty(peq), |
---|
329 | lhs = fun(1:(peq-1)); % left-hand side |
---|
330 | rhs = fun((peq+1):end); % right-hand side |
---|
331 | else |
---|
332 | lhs = ''; |
---|
333 | rhs = fun; |
---|
334 | end |
---|
335 | |
---|
336 | |
---|
337 | % process the lhs (if present) |
---|
338 | if ~isempty(lhs), |
---|
339 | pob = strfind(lhs,'('); % position of opening bracket |
---|
340 | pcb = strfind(lhs,')'); % position of closing bracket |
---|
341 | if ~isempty(pob) |
---|
342 | if pob==1, |
---|
343 | f.yvar = 'y'; |
---|
344 | else |
---|
345 | f.yvar = lhs(1:(pob-1)); |
---|
346 | end |
---|
347 | f.xvar = lhs((pob+1):(pcb-1)); |
---|
348 | else |
---|
349 | f.yvar = lhs; |
---|
350 | f.xvar = 'x'; |
---|
351 | end |
---|
352 | else % if no lhs present: |
---|
353 | f.xvar='x'; |
---|
354 | f.yvar='y'; |
---|
355 | end |
---|
356 | |
---|
357 | |
---|
358 | % process the 'poly' (polynomial fit) in the rhs |
---|
359 | if strfind(rhs,'poly'), % polynomial fit: |
---|
360 | order=str2double(rhs(5:end)); |
---|
361 | if isempty(order), % new v1.22 |
---|
362 | str_ord=inputdlg('Order of the polynomial fit','Polynomial order',1,{'2'}); |
---|
363 | if ~isempty(str_ord), |
---|
364 | order=str2double(str_ord{1}); |
---|
365 | f.name=['poly' str_ord{1}]; |
---|
366 | else |
---|
367 | clear f; |
---|
368 | return; |
---|
369 | end |
---|
370 | end |
---|
371 | if order>20, |
---|
372 | error('Ezyfit:ezfit:invalidPolynomDegree','Invalid polynom degree.'); |
---|
373 | end |
---|
374 | rhs = [fp.polynom_coeffname '0']; |
---|
375 | for i=1:order, |
---|
376 | if i>1 |
---|
377 | rhs = [rhs '+' fp.polynom_coeffname num2str(i) '*' f.xvar '^' num2str(i)]; |
---|
378 | else |
---|
379 | rhs = [rhs '+' fp.polynom_coeffname num2str(i) '*' f.xvar]; % new v2.22 |
---|
380 | end |
---|
381 | end |
---|
382 | end |
---|
383 | |
---|
384 | |
---|
385 | % search for option 'lin' or 'log': |
---|
386 | % (if several are present, use the last one) |
---|
387 | % (if none is present, check the Y-scale of the current figure) |
---|
388 | % (if no figure present, take 'lin'). |
---|
389 | fun = strrep([rhs ';' remfun], ';;', ';'); |
---|
390 | f.fitmode=''; |
---|
391 | plin=strfind(fun,';lin'); |
---|
392 | if ~isempty(plin), plin=plin(end); else plin=1; end |
---|
393 | plog=strfind(fun,';log'); |
---|
394 | if ~isempty(plog), plog=plog(end); else plog=1; end |
---|
395 | plast=max(plin,plog); |
---|
396 | if plast==1, % if nothing specified |
---|
397 | f.fitmode='lin'; % use 'lin' |
---|
398 | else |
---|
399 | f.fitmode=fun((plast+1):(plast+3)); |
---|
400 | end |
---|
401 | fun=strrep(fun,';lin',''); |
---|
402 | fun=strrep(fun,';log',''); |
---|
403 | |
---|
404 | |
---|
405 | % check for negative data in log mode (new v1.23, fixed 1.24): |
---|
406 | if (strcmp(f.fitmode,'log') && sum(y<=0)>0) |
---|
407 | disp('Warning: Zero or negative data ignored'); |
---|
408 | nonzero=find(y>0); |
---|
409 | x=x(nonzero); |
---|
410 | y=y(nonzero); |
---|
411 | end |
---|
412 | |
---|
413 | |
---|
414 | % separates again the first part (fitting function itself) and the |
---|
415 | % remainder: |
---|
416 | p = strfind(fun,';'); p=p(1); |
---|
417 | f.eq = fun(1:(p-1)); |
---|
418 | remfun = fun((p+1):end); |
---|
419 | |
---|
420 | |
---|
421 | % convert the equation in the matlab syntax (parameters named m(1), |
---|
422 | % m(2)...) |
---|
423 | [eqml,param] = eq2ml(f.eq, f.xvar); |
---|
424 | maxm = length(param); % number of parameters |
---|
425 | if maxm==0 % new v2.20 |
---|
426 | error('Ezyfit:ezfit:noParameter','No parameter to be fitted.'); |
---|
427 | end |
---|
428 | |
---|
429 | |
---|
430 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
431 | % processing the initial guess |
---|
432 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
433 | |
---|
434 | % initial guess for the m(i) by default (all m(i)=1): |
---|
435 | if ~exist('m0','var') |
---|
436 | m0=ones(1,maxm); |
---|
437 | switch usepredefinedfit % new v2.50 : "magic" initial guesses |
---|
438 | case 'linear' % m*x |
---|
439 | m0(1) = (y(end)-y(1))/(x(end)-x(1)); |
---|
440 | case 'affine' % a*x+b |
---|
441 | m0(1) = (y(end)-y(1))/(x(end)-x(1)); % a |
---|
442 | m0(2) = mean(y-m0(1)*x); % b |
---|
443 | case 'affineshift' % a*(x-b) |
---|
444 | m0(1) = (y(end)-y(1))/(x(end)-x(1)); % a |
---|
445 | m0(2) = mean(x-y/m0(1)); % b |
---|
446 | case 'power' % a*x^n; |
---|
447 | % if both x and y are of constant sign: |
---|
448 | if sum(diff(sign(x)))==0 || sum(diff(sign(y)))==0 |
---|
449 | m0(2) = log(y(end)/y(1)) / log(x(end)/x(1)); % n |
---|
450 | m0(1) = mean(y./(x.^m0(2))); % a |
---|
451 | end |
---|
452 | case 'powerc' % a*x^n+c; |
---|
453 | % if both x and y are of constant sign: |
---|
454 | if sum(diff(sign(x)))==0 || sum(diff(sign(y)))==0 |
---|
455 | m0(3) = log(y(end)/y(1)) / log(x(end)/x(1)); % n |
---|
456 | m0(1) = mean(y./(x.^m0(3))); % a |
---|
457 | end |
---|
458 | m0(2) = mean(y); % c |
---|
459 | case 'powershift' % a*(x+b)^n; |
---|
460 | % if both x and y are of constant sign: |
---|
461 | m0(2) = mean(x); % b |
---|
462 | m0(3) = log(y(end)/y(1)) / log((x(end)+m0(2))/(x(1)+m0(2))); % n |
---|
463 | m0(1) = mean(y./((x+m0(2)).^m0(3))); % a |
---|
464 | case 'exp' % a*exp(b*x) |
---|
465 | m0(2) = (log(y(end)/y(1))) / (x(end)-x(1)); % b |
---|
466 | m0(1) = mean(y./exp(m0(2)*x)); % a |
---|
467 | case 'expdiv' % a*exp(x/b) |
---|
468 | m0(2) = (x(end)-x(1)) / (log(y(end)/y(1))); % b |
---|
469 | m0(1) = mean(y./exp(x/m0(2))); % a |
---|
470 | case 'explim' % a*(1-exp(-x/b)) |
---|
471 | m0(1) = max(y); % a |
---|
472 | m0(2) = mean(-x./(log(1-y/m0(1)))); % b |
---|
473 | case 'expc' % a*exp(b*x)+c |
---|
474 | m0(3) = mean(y); % c |
---|
475 | m0(2) = (log((y(end)-m0(3))/(y(1)-m0(3)))) / (x(end)-x(1)); % b |
---|
476 | m0(1) = mean((y-m0(3))./exp(m0(2)*x)); % a |
---|
477 | case 'log' % a*log(b*x) |
---|
478 | m0(1) = (y(end)-y(1))/(log(x(end)/x(1))); % a |
---|
479 | m0(2) = mean(exp(y/m0(1))./x); % b |
---|
480 | case 'logc' % a*log(x)+b |
---|
481 | m0(1) = (y(end)-y(1))/(log(x(end)/x(1))); % a |
---|
482 | m0(2) = mean(y - m0(1)*log(x)); % b |
---|
483 | case {'sin','cos'} % a*sin(b*x), a*cos(b*x) |
---|
484 | m0(1) = std(y,1)*sqrt(2); % a |
---|
485 | m0(2) = 50/(x(end)-x(1)); % b |
---|
486 | case {'sinc','cosc'} % a*sin(b*x)+c, a*cos(b*x)+c |
---|
487 | m0(1) = std(y,1)*sqrt(2); % a |
---|
488 | m0(2) = 50/(x(end)-x(1)); % b |
---|
489 | m0(3) = mean(y); % c |
---|
490 | case 'sinphi' % a*sin(b*x+phi) |
---|
491 | m0(1) = std(y,1)*sqrt(2); % a |
---|
492 | m0(2) = 50/(x(end)-x(1)); % b |
---|
493 | m0(3) = 1; % phi |
---|
494 | case 'sinphic' % a*sin(b*x+phi)+c |
---|
495 | m0(1) = std(y,1)*sqrt(2); % a |
---|
496 | m0(2) = 50/(x(end)-x(1)); % b |
---|
497 | m0(3) = mean(y); % c |
---|
498 | m0(4) = 1; % phi |
---|
499 | case 'cngauss' |
---|
500 | m0(1) = (mean(x.^2.*y)/mean(y))^(1/2); % sigma |
---|
501 | case 'cfgauss' % a*exp(-(x^2)/(2*sigma^2)) |
---|
502 | m0(1) = max(y); % a |
---|
503 | m0(2) = (mean(x.^2.*y)/mean(y))^(1/2); % sigma |
---|
504 | case 'ngauss' % exp(-((x-x_0)^2)/(2*sigma^2))/(2*pi*sigma^2)^(1/2) |
---|
505 | m0(2) = mean(x.*y)/mean(y); % x_0 |
---|
506 | m0(1) = (mean((x-m0(2)).^2.*y)/mean(y))^(1/2); % sigma |
---|
507 | case {'fgauss','gauss'} % a*exp(-((x-x_0)^2)/(2*sigma^2)) |
---|
508 | m0(1) = max(y); % a |
---|
509 | m0(3) = mean(x.*y)/mean(y); % x_0 |
---|
510 | m0(2) = (mean((x-m0(3)).^2.*y)/mean(y))^(1/2); % sigma |
---|
511 | end |
---|
512 | % if an initial guess is 0 or infinity or has imaginary part, set it to 1 |
---|
513 | for np=1:length(m0) |
---|
514 | if m0(np)==0 || isinf(m0(np)) || imag(m0(np))~=0 |
---|
515 | m0(np) = 1; |
---|
516 | end |
---|
517 | end |
---|
518 | end |
---|
519 | |
---|
520 | |
---|
521 | % search for initial guess defined into remfun |
---|
522 | remfun = strrep([remfun ';'], ';;', ';'); % adds a final ';' |
---|
523 | while strfind(remfun,'='), |
---|
524 | peq=strfind(remfun,'='); peq=peq(1); |
---|
525 | pc=strfind(remfun,';'); pc=pc(1); |
---|
526 | if pc>peq, % if ';' after '=' |
---|
527 | par=remfun(1:(peq-1)); |
---|
528 | for i=1:maxm, |
---|
529 | if strcmp(param{i},par), |
---|
530 | m0(i)=eval(remfun((peq+1):(pc-1))); |
---|
531 | end |
---|
532 | end |
---|
533 | else |
---|
534 | error('Ezyfit:ezfit:syntaxError','Invalid syntax'); |
---|
535 | end |
---|
536 | remfun=remfun((pc+1):end); % removes the processed i.g. and loops back |
---|
537 | end |
---|
538 | |
---|
539 | |
---|
540 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
541 | % fitting itself |
---|
542 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
543 | |
---|
544 | switch f.fitmode, |
---|
545 | case 'lin', |
---|
546 | x__ref = x; |
---|
547 | try |
---|
548 | m=fminsearch(@fitlin, m0); |
---|
549 | catch |
---|
550 | error('Ezyfit:ezfit:fminsearchError','Fit: error during the fminsearch procedure'); |
---|
551 | end |
---|
552 | y_fit = eval(eqml); |
---|
553 | ssr = sum(abs(y_fit-mean(y)).^2); |
---|
554 | sse = sum(abs(y_fit-y).^2); |
---|
555 | f.r=sqrt(ssr/(sse+ssr)); |
---|
556 | case 'log', |
---|
557 | x__ref = x; |
---|
558 | try |
---|
559 | m=fminsearch(@fitlog, m0); % bug fixed here! (v2.51) |
---|
560 | catch |
---|
561 | error('Ezyfit:ezfit:fminsearchError','Fit: error during the fminsearch procedure'); |
---|
562 | end |
---|
563 | y_fit = eval(eqml); |
---|
564 | ssr = sum(abs(log(y_fit)-mean(log(y))).^2); |
---|
565 | sse = sum(abs(log(y_fit)-log(y)).^2); |
---|
566 | f.r=sqrt(ssr/(sse+ssr)); |
---|
567 | otherwise |
---|
568 | error('Ezyfit:ezfit:unknownMode''Unknown fit mode'); |
---|
569 | end |
---|
570 | |
---|
571 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
572 | % outputs |
---|
573 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
574 | |
---|
575 | |
---|
576 | % fills the output structure: |
---|
577 | f.param = param; |
---|
578 | f.m = m; |
---|
579 | f.m0 = m0; |
---|
580 | f.x = x; |
---|
581 | f.y = y; |
---|
582 | if sum(dy-1) % store the error bars only if defined |
---|
583 | f.dy = dy; |
---|
584 | end |
---|
585 | if exist('h','var'), f.hdata=h; end % handle to the data |
---|
586 | |
---|
587 | % stores the fit in the variable 'lastfit' in the 'base' workspace: |
---|
588 | assignin('base','lastfit',f); |
---|
589 | |
---|
590 | if strcmp(fp.automakevarfit,'on') |
---|
591 | makevarfit; |
---|
592 | end |
---|
593 | |
---|
594 | % ending displays (if no output argument): |
---|
595 | if ~nargout |
---|
596 | if strcmp(fp.dispeqmode,'on') % new v2.30 |
---|
597 | dispeqfit(f,fp); |
---|
598 | end |
---|
599 | clear f; |
---|
600 | end |
---|
601 | |
---|
602 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
603 | % End of the main function EZFIT |
---|
604 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
605 | |
---|
606 | |
---|
607 | |
---|
608 | % Nested functions that evaluate the fit for prescribed parameters m(i), |
---|
609 | % and return the chi2 (sum of the squared difference between the input |
---|
610 | % curve and the fit), in lin or log mode: |
---|
611 | |
---|
612 | function chi2 = fitlin(m) |
---|
613 | y_fit = eval(eqml); |
---|
614 | chi2 = sum(((y_fit - y).^2)./(dy.^2)); |
---|
615 | end |
---|
616 | |
---|
617 | function chi2 = fitlog(m) |
---|
618 | y_fit = eval(eqml); |
---|
619 | chi2 = sum((log(y_fit)-log(y)).^2); |
---|
620 | end |
---|
621 | end |
---|