www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregusermod/mmf.m

    function varargout= mmf(U,x,varargin)
%MMF Morgan-Mercer-Flodin user defined model for MBC
%
% varargout= mmf(U,x,varargin)
%   
% To use this function, the command 
%   U2= checkin(xregusermod,'funcname',xtest)
% must be run. The last argument is a column based input matrix to
% test the function. 

%  Copyright 2000-2014 The MathWorks, Inc. and Ford Global Technologies, Inc.



b= double(U);
if isa(x,'double')
    % shortcut for fast eval 
    
    % function definition
    % make sure this is vectorised
    y= b(1) - (b(1)-b(2))./(1+(b(3)*x).^b(4));
    y(x<0)=NaN;
    % this line must not be changed
    varargout{1}= y;
else
    %--------------------------------------------------------------------------
    % DO NOT CHANGE THIS CODE
    % x specifies which internal function to evaluate.
    % All internal functions are prefixed by 'i_'.
    % The model and model parameters are passed to all internal functions
    fInternal = str2func(sprintf('i_%s',lower(x)));
    [varargout{1:nargout}]= fInternal(U,b,varargin{:});
    %--------------------------------------------------------------------------
end

%---------------------------------------------------------------
% i_nfactors number of input factors 
%---------------------------------------------------------------
function n= i_nfactors(U,b) %#ok<INUSD,DEFNU>
% you must specify this 
n= 1;

%---------------------------------------------------------------
% i_numparams number of parameters 
%---------------------------------------------------------------
function n= i_numparams(U,b) %#ok<INUSD,DEFNU>
% you must specify this or initial
n= 4;

%---------------------------------------------------------------
% i_initial initial parameter values
%---------------------------------------------------------------
function [param,OK]= i_initial(U,p,X,Y) %#ok<INUSL,DEFNU>
% initial parameter values

if nargin==4
    % data dependent initial parameter values
    OK=1;
    % Bounds from spec.
    minbound=[2*eps eps eps eps]';
    maxbound=[inf inf inf inf]';
    
    if nargin==1 || isempty(X) || isempty(Y);
        param=[2 1 1 1]';  %default set
        return
    end
    Y=double(Y);
    X=double(X);
    
    % Remove any X,Y ==0
    xind= X==0;
    yind= Y==0;
    badind= or(xind,yind);
    Y(badind)=[];
    X(badind,:)=[];
    
    % Initial estimate of alpha
    mY=max(Y);
    r= mY-min(Y);
    % find out if there a number of points very close to the maximum y
    % we are then likely to be very close to the asymtote
    f=    min(length(find(Y-mY > -r*0.01))+1,5);
    alpha= (1 + 10^(-f))*mY;
    
    mY=min(Y);
    % find out if there a number of points very close to the maximum y
    % we are then likely to be very close to the asymtote
    f=    min(length(find(Y-mY < r*0.01))+1,5);
    beta= (1 - 10^(-f))*mY;
    
    % Create regression matrix
    Xs=[ones(size(X(:))) log(X(:))];
    
    % QR Decomposition
    %[Q,R]=qr(Xs,0);
    % Transform Y
    % Find residuals
    ylin= log((Y-beta)./(alpha-Y));
    %p= R\(Q'*ylin);
    p= Xs\ylin;
    
    delta=p(2);
    kappa=exp(p(1)/delta);
    
    
    % Returned parameter initial values
    param= [alpha beta kappa delta]';
    
    toosmall=find(param<minbound);
    toobig=find(param>maxbound);
    
    param(toosmall)=minbound(toosmall);
    param(toobig)=maxbound(toobig);
    
    OK= OK & isreal(param);
    
else
    param= [2 1 1 1]';
    OK=1;
end


% optional functions
% case 'jacobian',          % analytic jacobian
% case 'constraints'        % define parameter constraints
% case 'nlconstraints'      % evaluate nonlinear constraints
% case 'char'               % display string for function
% case 'str_func'           % one line summary of function
% case 'initial'            % initial parameter values initial(U,x,y)

% local regression optional functions

% case 'rfvals'             % response feature evaluation
% case 'reconstruct'        % local model reconstruction

%---------------------------------------------------------------
% i_constraints constraint definition for least squares fitting
%---------------------------------------------------------------
function [LB,UB,A,c,nlcon,optparams]= i_constraints(U,b,varargin) %#ok<INUSD,DEFNU>

% Lower Bounds, Upper Bounds  
% make sure these are column vectors
LB=[2*eps eps eps eps]';
UB=[inf inf inf inf]';

% only bounds are used for local regression at present

% Linear A, Linear b  (A*b<c)
A= [-1 1 0 0];
c= 0;
% Number of NonLinear constraints
nlcon= 0;

% Note that fmincon will be used if A and b is empty or nlcon>0
% otherwise lsqnonlin is used

% optional parameters for cost function
optparams= [];

%---------------------------------------------------------------
% i_foptions fitting options
%---------------------------------------------------------------
function fopts= i_foptions(U,b,fopts) %#ok<INUSL,DEFNU>

fopts= optimset(fopts,...
   'Display','none');

%---------------------------------------------------------------
% i_jacobian analytic jacobian (if defined)
%---------------------------------------------------------------
function J= i_jacobian(U,b,x) %#ok<INUSL,DEFNU>
% analytic jacobian (if defined)
x = x(:);

a=b(1);
bb=b(2);
k=b(3);
d=b(4);

kxd= (k*x).^d;

Ja= 1-1./(1+kxd);
Jb=1./(1+kxd);
Jk=(a-bb)./(1+kxd).^2.*kxd*d/k;
Jd=(a-bb)./(1+kxd).^2.*kxd.*log(k*x);

J= [Ja Jb Jk Jd];

% this line must not be changed
%---------------------------------------------------------------
% i_labels user parameter names 
%---------------------------------------------------------------
function c= i_labels(U,b) %#ok<INUSD,DEFNU>
% user parameter names 

c={'\alpha','\beta','\kappa','\delta'};

%---------------------------------------------------------------
% i_char display equation string
%---------------------------------------------------------------
function str= i_char(U,b) %#ok<INUSD,DEFNU>
% display equation string

s= get(U,'symbol');
% this can contain TeX expressions supported by HG text objects
str=sprintf('%.3g - (%.3g-%.3g)/(1+(%.3g*x)^{%.3g})',b([1 1 2 3 4]));
str = strrep(str,'/x^',['/',detex(s{1}),'^']);

%---------------------------------------------------------------
% i_str_func one line summary of function
%---------------------------------------------------------------
function str= i_str_func(U,b,TeX) %#ok<INUSL,DEFNU>
% one line summary of function
s= get(U,'symbol');
if nargin==2 || TeX
	s= detex(s);
end
lab= labels(U);

str = [lab{1},' - (',lab{1},'-',lab{2},')/(1+(',lab{3},'*',s{1},')^{',lab{4},'})'];

%---------------------------------------------------------------
% i_rfnames  response feature names 
%---------------------------------------------------------------
function [rname, default]= i_rfnames(U,b) %#ok<INUSD,DEFNU>
% response feature names 
rname= {'(Delta-1)/(2*Delta)'};

if nargout>1
    default = [1 2 3 5]; %{'alpha','beta','kappa','(Delta-1)/(2*Delta)'};
end

%---------------------------------------------------------------
% i_rfvals evaluate response features and gradient
%---------------------------------------------------------------
function [rf,dG]= i_rfvals(U,b) %#ok<INUSL,DEFNU>
% Note: model parameters are automatically available as response features

% this is an example of how to implement a nonlinear response feature
% response feature definition
rf= (b(4)-1)/(2*b(4));

if nargout>1
    % delrf/delbi
    dG=[0 0 0 1/(2*b(4)^2)];
end

%---------------------------------------------------------------
% i_reconstruct nonlinear reconstruction
%---------------------------------------------------------------
function p= i_reconstruct(U,b,Yrf,dG,rfuser) %#ok<INUSL,DEFNU>
% local reconstruction
% rfuser is an index to the user defined response features
% so you can figure out which rf's are which

% this solves for linear response features which can be estimated independently 
% of the nonlinear rf.
%  if all response features are linear you don't need to define 'reconstruct'
p= Yrf/dG';

% must have alpha, beta, kappa
% may have delta or the deltaRF
f=dG(:,4)~=0;
if ~any(rfuser==4)
    % we have the deltaRF
    p(:,end)= 1./(1-2*Yrf(:,f));
end

% k and delta must be +ve
p(:,3:4)= max(p(:,3:4),eps*16);

% this line must not be changed


%---------------------------------------------------------------
% i_evalbuild simulink evaluation block
%---------------------------------------------------------------
function Blk = i_evalbuild(U, b, sys) %#ok<INUSL,DEFNU>
% evalbuild method adds an evaluation simulink block to
% a simulink implementation of a model

Blk= add_block('mbcSLModels/LocalMod/mmfEval',[sys,'/mmfEval']);

% break library link
set_param(Blk,'linkstatus','none');


%---------------------------------------------------------------
% i_slrecon simulink reconstruction block
%---------------------------------------------------------------
function blk= i_slrecon(U, b, LUM, sname) %#ok<INUSL,DEFNU>
% GOMP/SLRECON -  adds the appropriate reconstruct block

blk= add_block('mbcSLModels/Reconstruct/mmfRecon',...
   [sname,'/Reconstruct']);

% break library link
set_param(blk,'linkstatus','none');