sensitivity copy.m

来自「MDPSAS工具箱是马里兰大学开发的」· M 代码 · 共 83 行

function [S, R] = sensitivity(A,Aparamkey);

% sensitivity.m Use finite-differences to compute sensitivity
%               of a naemodel object with respect to a parameter.
%
%               Called as
%
%               [S, R] = sensitivity(A,Aparamkey);

% Copyright (c) by Raymond A. Adomaitis, 1998-2005

% Based on the finite-difference algorithm for computing
% jacobian arrays (jacobian.m)

% check if residual in object A must be updated

if ~get(A,'residflag'), A = residual(A); end

if nargin ~= 2
   error('Must specify paramlist')
   end
   
paramindex = keyindex(get(A,'param'),Aparamkey);

if paramindex == 0, error('Invalid parameter'), end

R = columnnae(A,'resid');
P = get(A,'param',Aparamkey);

if isa(P,'scalarfield')
   P = get(P,'val');
   end

P = P(1:prod(size(P)),1);

Epsil = sqrt(eps)*max([abs(P),ones(size(P))],[],2).*sgn(P);
PPlusEpsil = P + Epsil;

% obsolete...
% Epsil = 1e-4*mean(abs(P));
% if Epsil < sqrt(eps), Epsil = sqrt(eps); end
% PPlusEpsil = P*(1+Epsil);
% PDiff = Epsil*P;
% I = find( abs(PDiff) < Epsil );
% PPlusEpsil(I) = P(I) + Epsil;
% PDiff(I) = Epsil;

n = 0;
B = A;

% Compute sensitivities of residuals to all keys in
% cell array Akey; each (set of) column(s) corresponds
% to each parameter in Akey.

for j = 1:length(P)

   B = set(B,'param',values( setval(get(A,'param'),PPlusEpsil(j),Aparamkey,j) ));
   B = residual(B);

   RPlusEpsil = columnnae(B,'resid');
   n = n+1;
   S(:,n) = ( RPlusEpsil - R )/Epsil(j);

% obsolete S(:,n) = ( RPlusEpsil - R )/PDiff(paramindex);

   end

% -------------------- sgn.m -------------------- %

function B = sgn(A)

% return the sign of A if nonzero; return 1 if zero

I = find( A==0 );
if ~isempty(I)
   B(I,1) = 1;
   end

J = find( A~=0 );
if ~isempty(J)
   B(J,1) = A(J,1)./abs(A(J,1));
   end