1function [nodes, weights] = cubature_with_gaussian_weight(d,n,method) % --*-- Unitary tests --*-- 2 3% Computes nodes and weights for a n-order cubature with gaussian weight. 4% 5% INPUTS 6% - d [integer] scalar, dimension of the region of integration. 7% - n [integer] scalar, approximation order (3 or 5). 8% - method [string] Method of approximation ('Stroud' or 'ScaledUnscentedTransform') 9% 10% OUTPUTS 11% - nodes [double] n×m matrix, with m=2×d if n=3 or m=2×d²+1 if n=5, nodes where the integrated function has to be evaluated. 12% - weights [double] m×1 vector, weights associated to the nodes. 13% 14% REMARKS 15% The routine returns nodes and associated weights to compute a multivariate integral of the form: 16% ∞ -<x,x> 17% ∫ f(x) × e dx 18% -∞ 19 20% Copyright (C) 2012-2019 Dynare Team 21% 22% This file is part of Dynare. 23% 24% Dynare is free software: you can redistribute it and/or modify 25% it under the terms of the GNU General Public License as published by 26% the Free Software Foundation, either version 3 of the License, or 27% (at your option) any later version. 28% 29% Dynare is distributed in the hope that it will be useful, 30% but WITHOUT ANY WARRANTY; without even the implied warranty of 31% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 32% GNU General Public License for more details. 33% 34% You should have received a copy of the GNU General Public License 35% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 36 37% Set default. 38if nargin<3 || isempty(method) 39 method = 'Stroud'; 40end 41 42if strcmp(method,'Stroud') && isequal(n,3) 43 r = sqrt(d); 44 nodes = r*[eye(d),-eye(d)]; 45 weights = ones(2*d,1)/(2*d); 46 return 47end 48 49if strcmp(method,'ScaledUnscentedTransform') && isequal(n,3) 50 % For alpha=1 and beta=kappa=0 we obtain the same weights and nodes than the 'Stroud' method (with n=3). 51 % For alpha=1, beta=0 and kappa=.5 we obtain sigma points with equal weights. 52 alpha = 1; 53 beta = 0; 54 kappa = 0.5; 55 lambda = (alpha^2)*(d+kappa) - d; 56 nodes = [ zeros(d,1) ( sqrt(d+lambda).*([ eye(d), -eye(d)]) ) ]; 57 w0_m = lambda/(d+lambda); 58 w0_c = w0_m + (1-alpha^2+beta); 59 weights = [w0_c; .5/(d+lambda)*ones(2*d,1)]; 60 return 61end 62 63if strcmp(method,'Stroud') && isequal(n,5) 64 r = sqrt((d+2)); 65 s = sqrt((d+2)/2); 66 m = 2*d^2+1; 67 A = 2/(n+2); 68 B = (4-d)/(2*(n+2)^2); 69 C = 1/(n+2)^2; 70 % Initialize the outputs 71 nodes = zeros(d,m); 72 weights = zeros(m,1); 73 % Set the weight for the first node (0) 74 weights(1) = A; 75 skip = 1; 76 % Set the remaining nodes and associated weights. 77 nodes(:,skip+(1:d)) = r*eye(d); 78 weights(skip+(1:d)) = B; 79 skip = skip+d; 80 nodes(:,skip+(1:d)) = -r*eye(d); 81 weights(skip+(1:d)) = B; 82 skip = skip+d; 83 for i=1:d-1 84 for j = i+1:d 85 nodes(:,skip+(1:4)) = s*ee(d,i,j); 86 weights(skip+(1:4)) = C; 87 skip = skip+4; 88 end 89 end 90 return 91end 92 93if strcmp(method,'Stroud') 94 error(['cubature_with_gaussian_weight:: Cubature (Stroud tables) is not yet implemented with n = ' int2str(n) '!']) 95end 96 97 98 99 100function v = e(n,i) 101v = zeros(n,1); 102v(i) = 1; 103 104function m = ee(n,i,j) 105m = zeros(n,4); 106m(:,1) = e(n,i)+e(n,j); 107m(:,2) = e(n,i)-e(n,j); 108m(:,3) = -m(:,2); 109m(:,4) = -m(:,1); 110 111return 112 113%@test:1 114d = 4; 115t = zeros(5,1); 116 117try 118 [nodes,weights] = cubature_with_gaussian_weight(d,3); 119 t(1) = 1; 120catch 121 t = t(1); 122 T = all(t); 123end 124 125if t(1) 126 m1 = nodes*weights; 127 m2 = nodes.^2*weights; 128 m3 = nodes.^3*weights; 129 m4 = nodes.^4*weights; 130 t(2) = dassert(m1,zeros(d,1),1e-12); 131 t(3) = dassert(m2,ones(d,1),1e-12); 132 t(4) = dassert(m3,zeros(d,1),1e-12); 133 t(5) = dassert(m4,d*ones(d,1),1e-10); 134 T = all(t); 135end 136%@eof:1 137 138%@test:2 139d = 4; 140Sigma = diag(1:d); 141Omega = diag(sqrt(1:d)); 142t = zeros(5,1); 143 144try 145 [nodes,weights] = cubature_with_gaussian_weight(d,3); 146 t(1) = 1; 147catch 148 t = t(1); 149 T = all(t); 150end 151 152if t(1) 153 nodes = Omega*nodes; 154 m1 = nodes*weights; 155 m2 = nodes.^2*weights; 156 m3 = nodes.^3*weights; 157 m4 = nodes.^4*weights; 158 t(2) = dassert(m1,zeros(d,1),1e-12); 159 t(3) = dassert(m2,transpose(1:d),1e-12); 160 t(4) = dassert(m3,zeros(d,1),1e-12); 161 t(5) = dassert(m4,d*transpose(1:d).^2,1e-10); 162 T = all(t); 163end 164%@eof:2 165 166%@test:3 167d = 4; 168Sigma = diag(1:d); 169Omega = diag(sqrt(1:d)); 170t = zeros(4,1); 171 172try 173 [nodes,weights] = cubature_with_gaussian_weight(d,3); 174 t(1) = 1; 175catch 176 t = t(1); 177 T = all(t); 178end 179 180if t(1) 181 nodes = Omega*nodes; 182 m1 = nodes*weights; 183 m2 = bsxfun(@times,nodes,transpose(weights))*transpose(nodes); 184 t(2) = dassert(m1,zeros(d,1),1e-12); 185 t(3) = dassert(diag(m2),transpose(1:d),1e-12); 186 t(4) = dassert(m2(:),vec(diag(diag(m2))),1e-12); 187 T = all(t); 188end 189%@eof:3 190 191%@test:4 192d = 10; 193a = randn(d,2*d); 194Sigma = a*a'; 195Omega = chol(Sigma,'lower'); 196t = zeros(4,1); 197 198try 199 [nodes,weights] = cubature_with_gaussian_weight(d,3); 200 t(1) = 1; 201catch 202 t = t(1); 203 T = all(t); 204end 205 206if t(1) 207 for i=1:length(weights) 208 nodes(:,i) = Omega*nodes(:,i); 209 end 210 m1 = nodes*weights; 211 m2 = bsxfun(@times,nodes,transpose(weights))*transpose(nodes); 212 m3 = nodes.^3*weights; 213 t(2) = dassert(m1,zeros(d,1),1e-12); 214 t(3) = dassert(m2(:),vec(Sigma),1e-12); 215 t(4) = dassert(m3,zeros(d,1),1e-12); 216 T = all(t); 217end 218%@eof:4 219 220%@test:5 221d = 5; 222t = zeros(6,1); 223 224try 225 [nodes,weights] = cubature_with_gaussian_weight(d,5); 226 t(1) = 1; 227catch 228 t = t(1); 229 T = all(t); 230end 231 232if t(1) 233 nodes = nodes; 234 m1 = nodes*weights; 235 m2 = nodes.^2*weights; 236 m3 = nodes.^3*weights; 237 m4 = nodes.^4*weights; 238 m5 = nodes.^5*weights; 239 t(2) = dassert(m1,zeros(d,1),1e-12); 240 t(3) = dassert(m2,ones(d,1),1e-12); 241 t(4) = dassert(m3,zeros(d,1),1e-12); 242 t(5) = dassert(m4,3*ones(d,1),1e-12); 243 t(6) = dassert(m5,zeros(d,1),1e-12); 244 T = all(t); 245end 246%@eof:5 247 248%@test:6 249d = 3; 250t = zeros(4,1); 251 252% Call the tested routine 253try 254 [nodes,weights] = cubature_with_gaussian_weight(d,3,'ScaledUnscentedTransform'); 255 t(1) = 1; 256catch 257 t = t(1); 258 T = all(t); 259end 260 261if t(1) 262 m1 = nodes*weights; 263 m2 = nodes.^2*weights; 264 m3 = nodes.^3*weights; 265 t(2) = dassert(m1,zeros(d,1),1e-12); 266 t(3) = dassert(m2,ones(d,1),1e-12); 267 t(4) = dassert(m3,zeros(d,1),1e-12); 268 T = all(t); 269end 270%@eof:6