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