Morris_Measure_Groups

 [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)

 Given the Morris sample matrix, the output values and the group matrix compute the Morris measures
 -------------------------------------------------------------------------
 INPUTS
 -------------------------------------------------------------------------
 Group [NumFactor, NumGroups] := Matrix describing the groups. 
 Each column represents one group. 
 The element of each column are zero if the factor is not in the
 group. Otherwise it is 1.

 Sample := Matrix of the Morris sampled trajectories 

 Output := Matrix of the output(s) values in correspondence of each point
 of each trajectory

 k = Number of factors
 -------------------------------------------------------------------------
 OUTPUTS 
 OutMatrix (NumFactor*NumOutputs, 3)= [Mu*, Mu, StDev]
 for each output it gives the three measures of each factor
 -------------------------------------------------------------------------

0001 function [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)
0002 % [SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group)
0003 %
0004 % Given the Morris sample matrix, the output values and the group matrix compute the Morris measures
0005 % -------------------------------------------------------------------------
0006 % INPUTS
0007 % -------------------------------------------------------------------------
0008 % Group [NumFactor, NumGroups] := Matrix describing the groups.
0009 % Each column represents one group.
0010 % The element of each column are zero if the factor is not in the
0011 % group. Otherwise it is 1.
0012 %
0013 % Sample := Matrix of the Morris sampled trajectories
0014 %
0015 % Output := Matrix of the output(s) values in correspondence of each point
0016 % of each trajectory
0017 %
0018 % k = Number of factors
0019 % -------------------------------------------------------------------------
0020 % OUTPUTS
0021 % OutMatrix (NumFactor*NumOutputs, 3)= [Mu*, Mu, StDev]
0022 % for each output it gives the three measures of each factor
0023 % -------------------------------------------------------------------------
0024 
0025 % Copyright (C) 2012 Dynare Team
0026 %
0027 % This file is part of Dynare.
0028 %
0029 % Dynare is free software: you can redistribute it and/or modify
0030 % it under the terms of the GNU General Public License as published by
0031 % the Free Software Foundation, either version 3 of the License, or
0032 % (at your option) any later version.
0033 %
0034 % Dynare is distributed in the hope that it will be useful,
0035 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0036 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0037 % GNU General Public License for more details.
0038 %
0039 % You should have received a copy of the GNU General Public License
0040 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0041 
0042 if nargin==0,
0043   disp(' ')
0044   disp('[SAmeas, OutMatrix] = Morris_Measure_Groups(NumFact, Sample, Output, p, Group);')
0045   return
0046 end
0047 
0048 OutMatrix=[];
0049 if nargin < 5, Group=[]; end
0050 
0051 NumGroups = size(Group,2);
0052 if nargin < 4 | isempty(p),
0053     p = 4;
0054 end
0055 Delt = p/(2*p-2);
0056 
0057 if NumGroups ~ 0
0058     sizea = NumGroups;      % Number of groups
0059     GroupMat=Group;
0060     GroupMat = GroupMat';
0061 else 
0062     sizea = NumFact; 
0063 end
0064 r=size(Sample,1)/(sizea+1);     % Number of trajectories
0065 
0066 % For Each Output
0067 for k=1:size(Output,2)
0068     
0069     OutValues=Output(:,k);
0070   
0071     % For each r trajectory
0072     for i=1:r
0073     
0074         % For each step j in the trajectory
0075         % Read the orientation matrix fact for the r-th sampling
0076         % Read the corresponding output values
0077         Single_Sample = Sample(i+(i-1)*sizea:i+(i-1)*sizea+sizea,:); 
0078         Single_OutValues = OutValues(i+(i-1)*sizea:i+(i-1)*sizea+sizea,:); 
0079         A = (Single_Sample(2:sizea+1,:)-Single_Sample(1:sizea,:))';
0080         Delta = A(find(A));
0081 
0082         % For each point of the fixed trajectory compute the values of the Morris function. The function
0083         % is partitioned in four parts, from order zero to order 4th.
0084         for j=1:sizea   % For each point in the trajectory i.e for each factor
0085             % matrix of factor which changes
0086             if NumGroups ~ 0;
0087                 AuxFind (:,1) = A(:,j);
0088 %                 AuxFind(find(A(:,j)),1)=1;
0089 %                 Pippo = sum((Group - repmat(AuxFind,1,NumGroups)),1);
0090 %                 Change_factor(j,i) = find(Pippo==0);
0091                 Change_factor = find(abs(AuxFind)>1e-010); 
0092                 % If we deal with groups we can only estimate the new mu*
0093                 % measure since factors in the same groups can move in
0094                 % opposite direction and the definition of the standard
0095                 % Morris mu cannopt be applied.
0096                 % In the new version the elementary effect is defined with
0097                 % the absolute value.
0098                 %SAmeas(find(GroupMat(Change_factor(j,i),:)),i) = abs((Single_OutValues(j) - Single_OutValues(j+1) )/Delt); %(2/3));
0099                 SAmeas(i,Change_factor') = abs((Single_OutValues(j) - Single_OutValues(j+1) )/Delt);   
0100             else
0101                 Change_factor(j,i) = find(Single_Sample(j+1,:)-Single_Sample(j,:));
0102                 % If no groups --> we compute both the original and
0103                 % modified measure
0104                 if Delta(j) > 0                              %=> +Delta
0105                     SAmeas(Change_factor(j,i),i) = (Single_OutValues(j+1) - Single_OutValues(j) )/Delt; %(2/3);
0106                 else                                         %=> -Delta
0107                     SAmeas(Change_factor(j,i),i) = (Single_OutValues(j) - Single_OutValues(j+1) )/Delt; %(2/3);
0108                 end 
0109             end
0110         end   %for j=1:sizea
0111     
0112     end     %for i=1:r
0113    
0114     if NumGroups ~ 0
0115         SAmeas = SAmeas';
0116     end
0117 
0118     % Compute Mu AbsMu and StDev
0119     if any(any(isnan(SAmeas)))
0120       for j=1:NumFact,
0121         SAm = SAmeas(j,:);
0122         SAm = SAm(find(~isnan(SAm)));
0123         rr=length(SAm);
0124         AbsMu(j,1) = sum(abs(SAm),2)/rr;
0125       if NumGroups == 0
0126         Mu(j,1) = sum(SAm,2)/rr;
0127         StDev(j,1) = sum((SAm - repmat(Mu(j),1,rr)).^2/(rr*(rr-1)),2).^0.5;
0128       end
0129       end
0130     else
0131       AbsMu = sum(abs(SAmeas),2)/r;
0132       if NumGroups == 0
0133         Mu = sum(SAmeas,2)/r;
0134         StDev = sum((SAmeas - repmat(Mu,1,r)).^2/(r*(r-1)),2).^0.5;
0135       end
0136     end
0137 
0138     % Define the output Matrix - if we have groups we cannot define the old
0139     % measure mu, only mu* makes sense
0140     if NumGroups > 0
0141         OutMatrix = [OutMatrix; AbsMu];   
0142     else
0143         OutMatrix = [OutMatrix; AbsMu, Mu, StDev];   
0144     end
0145 end     % For Each Output