[Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat)
Inputs: k (1,1) := number of factors examined or number of groups examined.
In case the groups are chosen the number of factors is stores in NumFact and
sizea becomes the number of created groups.
NumFact (1,1) := number of factors examined in the case when groups are chosen
r (1,1) := sample size
p (1,1) := number of intervals considered in [0, 1]
UB(sizea,1) := Upper Bound for each factor
LB(sizea,1) := Lower Bound for each factor
GroupNumber(1,1) := Number of groups (eventually 0)
GroupMat(NumFact,GroupNumber):= Matrix which describes the chosen groups. Each column represents a group and its elements
are set to 1 in correspondence of the factors that belong to the fixed group. All
the other elements are zero.
Local Variables:
sizeb (1,1) := sizea+1
sizec (1,1) := 1
randmult (sizea,1) := vector of random +1 and -1
perm_e(1,sizea) := vector of sizea random permutated indeces
fact(sizea) := vector containing the factor varied within each traj
DDo(sizea,sizea) := D* in Morris, 1991
A(sizeb,sizea) := Jk+1,k in Morris, 1991
B(sizeb,sizea) := B in Morris, 1991
Po(sizea,sizea) := P* in Morris, 1991
Bo(sizeb,sizea) := B* in Morris, 1991
Ao(sizeb,sizec) := Jk+1,1 in Morris, 1991
xo(sizec,sizea) := x* in Morris, 1991 (starting point for the trajectory)
In(sizeb,sizea) := for each loop orientation matrix. It corresponds to a trajectory
of k step in the parameter space and it provides a single elementary
effect per factor
MyInt()
Fact(sizea,1) := for each loop vector indicating which factor or group of factors has been changed
in each step of the trajectory
AuxMat(sizeb,sizea) := Delta*0.5*((2*B - A) * DD0 + A) in Morris, 1991. The AuxMat is used as in Morris design
for single factor analysis, while it constitutes an intermediate step for the group analysis.
Output: Outmatrix(sizeb*r, sizea) := for the entire sample size computed In(i,j) matrices
OutFact(sizea*r,1) := for the entire sample size computed Fact(i,1) vectors
Note: B0 is constructed as in Morris design when groups are not considered. When groups are considered the routine
follows the following steps:
1- Creation of P0 and DD0 matrices defined in Morris for the groups. This means that the dimensions of these
2 matrices are (GroupNumber,GroupNumber).
2- Creation of AuxMat matrix with (GroupNumber+1,GroupNumber) elements.
3- Definition of GroupB0 starting from AuxMat, GroupMat and P0.
4- The final B0 for groups is obtained as [ones(sizeb,1)*x0' + GroupB0]. The P0 permutation is present in GroupB0
and it's not necessary to permute the matrix (ones(sizeb,1)*x0') because it's already randomly created.
Reference:
A. Saltelli, K. Chan, E.M. Scott, "Sensitivity Analysis" on page 68 ss
F. Campolongo, J. Cariboni, JRC - IPSC Ispra, Varese, IT
Last Update: 15 November 2005 by J.Cariboni0001 function [Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat) 0002 %[Outmatrix, OutFact] = Sampling_Function_2(p, k, r, UB, LB, GroupMat) 0003 % Inputs: k (1,1) := number of factors examined or number of groups examined. 0004 % In case the groups are chosen the number of factors is stores in NumFact and 0005 % sizea becomes the number of created groups. 0006 % NumFact (1,1) := number of factors examined in the case when groups are chosen 0007 % r (1,1) := sample size 0008 % p (1,1) := number of intervals considered in [0, 1] 0009 % UB(sizea,1) := Upper Bound for each factor 0010 % LB(sizea,1) := Lower Bound for each factor 0011 % GroupNumber(1,1) := Number of groups (eventually 0) 0012 % GroupMat(NumFact,GroupNumber):= Matrix which describes the chosen groups. Each column represents a group and its elements 0013 % are set to 1 in correspondence of the factors that belong to the fixed group. All 0014 % the other elements are zero. 0015 % Local Variables: 0016 % sizeb (1,1) := sizea+1 0017 % sizec (1,1) := 1 0018 % randmult (sizea,1) := vector of random +1 and -1 0019 % perm_e(1,sizea) := vector of sizea random permutated indeces 0020 % fact(sizea) := vector containing the factor varied within each traj 0021 % DDo(sizea,sizea) := D* in Morris, 1991 0022 % A(sizeb,sizea) := Jk+1,k in Morris, 1991 0023 % B(sizeb,sizea) := B in Morris, 1991 0024 % Po(sizea,sizea) := P* in Morris, 1991 0025 % Bo(sizeb,sizea) := B* in Morris, 1991 0026 % Ao(sizeb,sizec) := Jk+1,1 in Morris, 1991 0027 % xo(sizec,sizea) := x* in Morris, 1991 (starting point for the trajectory) 0028 % In(sizeb,sizea) := for each loop orientation matrix. It corresponds to a trajectory 0029 % of k step in the parameter space and it provides a single elementary 0030 % effect per factor 0031 % MyInt() 0032 % Fact(sizea,1) := for each loop vector indicating which factor or group of factors has been changed 0033 % in each step of the trajectory 0034 % AuxMat(sizeb,sizea) := Delta*0.5*((2*B - A) * DD0 + A) in Morris, 1991. The AuxMat is used as in Morris design 0035 % for single factor analysis, while it constitutes an intermediate step for the group analysis. 0036 % 0037 % Output: Outmatrix(sizeb*r, sizea) := for the entire sample size computed In(i,j) matrices 0038 % OutFact(sizea*r,1) := for the entire sample size computed Fact(i,1) vectors 0039 % 0040 % Note: B0 is constructed as in Morris design when groups are not considered. When groups are considered the routine 0041 % follows the following steps: 0042 % 1- Creation of P0 and DD0 matrices defined in Morris for the groups. This means that the dimensions of these 0043 % 2 matrices are (GroupNumber,GroupNumber). 0044 % 2- Creation of AuxMat matrix with (GroupNumber+1,GroupNumber) elements. 0045 % 3- Definition of GroupB0 starting from AuxMat, GroupMat and P0. 0046 % 4- The final B0 for groups is obtained as [ones(sizeb,1)*x0' + GroupB0]. The P0 permutation is present in GroupB0 0047 % and it's not necessary to permute the matrix (ones(sizeb,1)*x0') because it's already randomly created. 0048 % Reference: 0049 % A. Saltelli, K. Chan, E.M. Scott, "Sensitivity Analysis" on page 68 ss 0050 % 0051 % F. Campolongo, J. Cariboni, JRC - IPSC Ispra, Varese, IT 0052 % Last Update: 15 November 2005 by J.Cariboni 0053 0054 % Copyright (C) 2012 Dynare Team 0055 % 0056 % This file is part of Dynare. 0057 % 0058 % Dynare is free software: you can redistribute it and/or modify 0059 % it under the terms of the GNU General Public License as published by 0060 % the Free Software Foundation, either version 3 of the License, or 0061 % (at your option) any later version. 0062 % 0063 % Dynare is distributed in the hope that it will be useful, 0064 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0065 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0066 % GNU General Public License for more details. 0067 % 0068 % You should have received a copy of the GNU General Public License 0069 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0070 0071 % Parameters and initialisation of the output matrix 0072 sizea = k; 0073 Delta = p/(2*p-2); 0074 %Delta = 1/3 0075 NumFact = sizea; 0076 GroupNumber = size(GroupMat,2); 0077 0078 if GroupNumber ~ 0; 0079 sizea = size(GroupMat,2); 0080 end 0081 0082 sizeb = sizea + 1; 0083 sizec = 1; 0084 Outmatrix = []; 0085 OutFact = []; 0086 0087 % For each i generate a trajectory 0088 for i=1:r 0089 0090 % Construct DD0 - OLD VERSION - it does not need communication toolbox 0091 % RAND(N,M) is an NXM matrix with random entries, chosen from a uniform distribution on the interval (0.0,1.0). 0092 % Note that DD0 tells if the factor have to be increased or ddecreased 0093 % by Delta. 0094 randmult = ones(k,1); 0095 v = rand(k,1); 0096 randmult (find(v < 0.5))=-1; 0097 randmult = repmat(randmult,1,k); 0098 DD0 = randmult .* eye(k); 0099 0100 % Construct DD0 - NEW VERSION - it needs communication toolbox 0101 % randsrc(m) generates an m-by-m matrix, each of whose entries independently takes the value -1 with probability 1/2, 0102 % and 1 with probability 1/2. 0103 % DD0 = randsrc(NumFact) .* eye(NumFact); 0104 0105 % Construct B (lower triangular) 0106 B = ones(sizeb,sizea); 0107 for j = 1:sizea 0108 B(1:j,j)=0; 0109 end 0110 0111 % Construct A0, A 0112 A0 = ones(sizeb,1); 0113 A = ones(sizeb,NumFact); 0114 0115 % Construct the permutation matrix P0. In each column of P0 one randomly chosen element equals 1 0116 % while all the others equal zero. 0117 % P0 tells the order in which order factors are changed in each 0118 % trajectory. P0 is created as it follows: 0119 % 1) All the elements of P0 are set equal to zero ==> P0 = zeros (sizea, sizea); 0120 % 2) The function randperm create a random permutation of integer 1:sizea, without repetitions ==> perm_e; 0121 % 3) In each column of P0 the element indicated in perm_e is set equal to one. 0122 % Note that P0 is then used reading it by rows. 0123 P0 = zeros (sizea, sizea); 0124 perm_e = randperm(sizea); % RANDPERM(n) is a random permutation of the integers from 1 to n. 0125 for j = 1:sizea 0126 P0(perm_e(j),j) = 1; 0127 end 0128 0129 % When groups are present the random permutation is done only on B. The effect is the same since 0130 % the added part (A0*x0') is completely random. 0131 if GroupNumber ~ 0 0132 B = B * (GroupMat*P0')'; 0133 end 0134 0135 % Compute AuxMat both for single factors and groups analysis. For Single factors analysis 0136 % AuxMat is added to (A0*X0) and then permutated through P0. When groups are active AuxMat is 0137 % used to build GroupB0. AuxMat is created considering DD0. If the element on DD0 diagonal 0138 % is 1 then AuxMat will start with zero and add Delta. If the element on DD0 diagonal is -1 0139 % then DD0 will start Delta and goes to zero. 0140 AuxMat = Delta*0.5*((2*B - A) * DD0 + A); 0141 0142 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0143 % a --> Define the random vector x0 for the factors. Note that x0 takes value in the hypercube 0144 % [0,...,1-Delta]*[0,...,1-Delta]*[0,...,1-Delta]*[0,...,1-Delta] 0145 MyInt = repmat([0:(1/(p-1)):(1-Delta)],NumFact,1); % Construct all possible values of the factors 0146 0147 % OLD VERSION - it needs communication toolbox 0148 % w = randint(NumFact,1,[1,size(MyInt,2)]); 0149 0150 % NEW VERSION - construct a version of random integers 0151 % 1) create a vector of random integers 0152 % 2) divide [0,1] into the needed steps 0153 % 3) check in which interval the random numbers fall 0154 % 4) generate the corresponding integer 0155 v = repmat(rand(NumFact,1),1,size(MyInt,2)+1); % 1) 0156 IntUsed = repmat([0:1/size(MyInt,2):1],NumFact,1); % 2) 0157 DiffAuxVec = IntUsed - v; % 3) 0158 0159 for ii = 1:size(DiffAuxVec,1) 0160 w(1,ii) = max(find(DiffAuxVec(ii,:)<0)); % 4) 0161 end 0162 x0 = MyInt(1,w)'; % Define x0 0163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0164 % b --> Compute the matrix B*, here indicated as B0. Each row in B0 is a 0165 % trajectory for Morris Calculations. The dimension of B0 is (Numfactors+1,Numfactors) 0166 if GroupNumber ~ 0 0167 B0 = (A0*x0' + AuxMat); 0168 else 0169 B0 = (A0*x0' + AuxMat)*P0; 0170 end 0171 0172 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0173 % c --> Compute values in the original intervals 0174 % B0 has values x(i,j) in [0, 1/(p -1), 2/(p -1), ... , 1]. 0175 % To obtain values in the original intervals [LB, UB] we compute 0176 % LB(j) + x(i,j)*(UB(j)-LB(j)) 0177 In = repmat(LB,1,sizeb)' + B0 .* repmat((UB-LB),1,sizeb)'; 0178 0179 % Create the Factor vector. Each component of this vector indicate which factor or group of factor 0180 % has been changed in each step of the trajectory. 0181 for j=1:sizea 0182 Fact(1,j) = find(P0(j,:)); 0183 end 0184 Fact(1,sizea+1) = 0; 0185 0186 Outmatrix = [Outmatrix; In]; 0187 OutFact = [OutFact; Fact']; 0188 0189 end