曲线拟合

nshukwrd

1. function [x,OPTIONS,CostFunction,JACOB] = curvefit(FUN,x0,XDATA,YDATA,OPTIONS,GRADFUN,varargin)
   
2. %曲线拟合
   
3. %   C=CURVEFIT('拟合模型',参数C初值,X数据,Y数据)
   
4. %求参数C使得 sum {(FUN(C,X数据)-Y数据).^2}最小化
   
5. %例
   
6. %   x | 0.1  0.2  0.15 0.0  -0.2 0.3
   
7. %   --|------------------------------
   
8. %   y | 0.95 0.84 0.86 1.06 1.50 0.72
   
9. %拟合 y=a*exp(b*x).先写M-函数fitfun2.m
   
10. %                   function   f=fitfun2(c,x)
    
11. %                   f=c(1)*exp(c(2)*x);
    
12. %  解法
    
13. %  x=[0.1,0.2,0.15,0,-0.2,0.3];
    
14. %  y=[0.95,0.84,0.86,1.06,1.50,0.72];
    
15. %  c=curvefit('fitfun2',[1,1]',x,y);
    
16. %  a=c(1),b=c(2)
    
17. %  f='a*exp(b*xi)';
    
18. %  xi=-0.2:0.01:0.3;
    
19. %  yi=eval(f);
    
20. %  plot(x,y,'o',xi,yi,'k')
    
21. %  title('CURVEFIT');   
    
22. 
    
23. %CURVEFIT Solves non-linear least squares problems.
    
24. %   CURVEFIT solves problems of the form:
    
25. %   min  sum {(FUN(X,XDATA)-YDATA).^2}  where FUN, XDATA and YDATA are
    
26. %    X                                  matrices.   
    
27. %
    
28. %   X=CURVEFIT('FUN',X0,XDATA,YDATA) starts at X0 and finds
    
29. %   coefficients X to best fit the nonlinear function FUN(X,XDATA)
    
30. %   to the data YDATA (in the least-squares sense). FUN is an M-file
    
31. %   that computes a function of X and XDATA and returns a matrix of
    
32. %   the objective function values: F=FUN(X,XDATA).
    
33. %   NOTE: YDATA must be the same size as the matrix F returned by FUN.
    
34. %
    
35. %   X=CURVEFIT('FUN',X0,XDATA,YDATA,OPTIONS) allows a vector of optional
    
36. %   parameters to be defined. OPTIONS(2) is a measure of the precision
    
37. %   required for the values of X at the solution. OPTIONS(3) is a measure
    
38. %   of the precision required of the objective function at the solution.
    
39. %   See HELP FOPTIONS.
    
40. %
    
41. %   X=CURVEFIT('FUN',X0,XDATA,YDATA,OPTIONS,'GRADFUN') enables a
    
42. %   function 'GRADFUN' to be entered which returns the partial derivatives
    
43. %   of the functions, dF/dX, (stored in columns) at the point X:
    
44. %   gf = GRADFUN(X,XDATA).
    
45. %
    
46. %   X=CURVEFIT('FUN',X,XDATA,YDATA,OPTIONS,'GRADFUN',P1,P2,..) passes the
    
47. %   problem-dependent parameters P1,P2,... directly to the functions FUN
    
48. %   and GRADFUN: FUN(X,XDATA,P1,P2,...) and GRADFUN(X,XDATA,P1,P2,...).  
    
49. %   Pass empty matrices for OPTIONS and 'GRADFUN' to use the default values.
    
50. %
    
51. %   [X,OPTIONS,F,J]=CURVEFIT('FUN',X0,XDATA,YDATA,...) returns, F,
    
52. %   the value of FUN(X,XDATA)-YDATA at the solution X, and J the Jacobian
    
53. %   of the function FUN at the solution.
    
54. %
    
55. %   FUN must be an M-file and not an inline object or expression. Use LEASTSQ
    
56. %   instead on inline objects or expressions.
    
57. 
    
58. %   Copyright (c) 1990-98 by The MathWorks, Inc.
    
59. %   $Revision: 1.6 $  $Date: 1997/11/29 01:23:03 $
    
60. %   Mary Ann Branch 8-22-96.
    
61. 
    
62. %   The default algorithm is the Levenberg-Marquardt method with a
    
63. %   mixed quadratic and cubic line search procedure.  A Gauss-Newton
    
64. %   method is selected by setting  OPTIONS(5)=1.
    
65. %
    
66. 
    
67. % ------------Initialization----------------
    
68. 
    
69. if nargin < 4, error('curvefit requires four input arguments');end
    
70. if nargin < 5, OPTIONS=[]; end
    
71. if nargin < 6, GRADFUN=[]; end
    
72. 
    
73. % Need extra argname for XDATA
    
74. curvefitarg = 1;
    
75. lenVarIn = length(varargin);
    
76. 
    
77. % Convert to inline function as needed.
    
78. if isempty(FUN)
    
79.    error('FUN must be a function name.');
    
80. elseif isa(FUN,'inline')
    
81.    error('FUN must be a function name, not an inline object. Use LEASTSQ instead.');
    
82. elseif isstr(FUN)
    
83.    % Make sure FUN isn't an expression via fcnchk
    
84.    fun0 = fcnchk(FUN);
    
85.    if isa(fun0,'inline')
    
86.       error('FUN must be a function name, not a string expression. Use LEASTSQ instead.');
    
87.    else
    
88.       % Momentarily turn FUN into an inline to get at the argnames and formula
    
89.       fun1 = inline(FUN,lenVarIn+1);
    
90.       fun1args = argnames(fun1);
    
91.       % Compute the inline function that is the difference with YDATA      
    
92.       arglist = sprintf('%s(', formula(fun1));
    
93.       for k = 1:(nargin(fun1)-1);
    
94.          arglist = sprintf('%s%s,', arglist, deblank(fun1args{k,:}));
    
95.       end
    
96.       arglist = sprintf('%s%s)',arglist, deblank(fun1args{nargin(fun1),:}));
    
97.       funfcn = inline([arglist, ...
    
98.             '-P',int2str(lenVarIn+2)],lenVarIn+2);
    
99.    end
    
100. else
     
101.    error('FUN is an unrecognized data type.')
     
102. end
     
103. 
     
104. if ~isempty(GRADFUN)
     
105.    if isa(GRADFUN,'inline')
     
106.       error('GRADFUN must be a function name and not an inline object.');
     
107.    elseif isstr(GRADFUN)
     
108.       % Make sure FUN isn't an expression via fcnchk
     
109.       fun0 = fcnchk(GRADFUN);
     
110.       if isa(fun0,'inline')
     
111.          error('GRADFUN must be a function name, not a string expression.');
     
112.       else
     
113.          % Momentarily turn GRADFUN into an inline to get at the argnames and formula
     
114.          fun1 = inline(GRADFUN,lenVarIn+1);
     
115.          fun1args = argnames(fun1);
     
116.          % Compute the inline function that is the difference with YDATA
     
117.          arglist = sprintf('%s(', formula(fun1));
     
118.          for k = 1:(nargin(fun1)-1);
     
119.             arglist = sprintf('%s%s,', arglist, deblank(fun1args{k,:}));
     
120.          end
     
121.          arglist = sprintf('%s%s)',arglist, deblank(fun1args{nargin(fun1),:}));
     
122.          gradfcn = inline(arglist,lenVarIn+2);
     
123.       end   
     
124.    else
     
125.       error('GRADFUN is an unrecognized data type.')
     
126.    end
     
127. else
     
128.    gradfcn = [];
     
129. end
     
130. 
     
131. 
     
132. % Check the size of f and YDATA
     
133. try
     
134.    f = feval(FUN,x0,XDATA,varargin{:});
     
135. catch
     
136.    error(['Error in evaluating user supplied function: ',FUN])
     
137. end
     
138. 
     
139. if ~isequal(size(f), size(YDATA))
     
140.    error('Function value and YDATA sizes are incommensurate.')
     
141. end
     
142. 
     
143. % The order YDATA, XDATA is relevant: match P2 to YDATA
     
144. [x,OPTIONS,CostFunction,JACOB] = ...
     
145.    nlsqold(funfcn,x0,OPTIONS,gradfcn,XDATA,varargin{:},YDATA);
     
146. 
     
147. %--end of curvefit--
     
148. 
     

复制代码