1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
% STK_SAMPLING_SOBOL generates points from a Sobol sequence
%
% CALL: X = stk_sampling_sobol (N, D)
%
% computes the first N terms of a D-dimensional Sobol sequence (with
% N < 2^32 and D <= 1111). The sequence is generated using the algorithm
% of Bratley and Fox [1], as modified by Joe and Kuo [3].
%
% CALL: X = stk_sampling_sobol (N, DIM, BOX)
%
% does the same thing in the DIM-dimensional hyperrectangle specified by the
% argument BOX, which is a 2 x DIM matrix where BOX(1, j) and BOX(2, j) are
% the lower- and upper-bound of the interval on the j^th coordinate.
%
% CALL: X = stk_sampling_sobol (N, D, BOX, DO_SKIP)
%
% skips an initial segment of the Sobol sequence if DO_SKIP is true. More
% precisely, according to the recommendation of [2] and [3], a number of
% points equal to the largest power of 2 smaller than n is skipped. If
% DO_SKIP is false, the beginning of the sequence is returns, as in the
% previous cases (in other words, DO_SKIP = false is the default).
%
% NOTE: Implementation
%
% The C implementation under the hood is due to Steven G. Johnson, and
% was borrowed from the NLopt toolbox [4] (version 2.4.2).
%
% REFERENCE
%
% [1] Paul Bratley and Bennett L. Fox, "Algorithm 659: Implementing Sobol's
% quasirandom sequence generator", ACM Transactions on Mathematical
% Software, 14(1):88-100, 1988.
%
% [2] Peter Acworth, Mark Broadie and Paul Glasserman, "A Comparison of Some
% Monte Carlo and Quasi Monte Carlo Techniques for Option Pricing", in
% Monte Carlo and Quasi-Monte Carlo Methods 1996, Lecture Notes in
% Statistics 27:1-18, Springer, 1998.
%
% [3] Stephen Joe and Frances Y. Kuo, "Remark on Algorithm 659: Implementing
% Sobol's Quasirandom Sequence Generator', ACM Transactions on
% Mathematical Software, 29(1):49-57, 2003.
%
% [4] Steven G. Johnson, The NLopt nonlinear-optimization package,
% http://ab-initio.mit.edu/nlopt.
%
% SEE ALSO: stk_sampling_halton_rr2
% Copyright Notice
%
% Copyright (C) 2016-2018 CentraleSupelec
%
% Author: Julien Bect <julien.bect@centralesupelec.fr>
% Copying Permission Statement
%
% This file is part of
%
% STK: a Small (Matlab/Octave) Toolbox for Kriging
% (https://github.com/stk-kriging/stk/)
%
% STK is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your
% option) any later version.
%
% STK is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with STK. If not, see <http://www.gnu.org/licenses/>.
function x = stk_sampling_sobol (n, dim, box, do_skip)
% Default values
if nargin < 4
do_skip = false;
if nargin < 3
box = [];
if nargin < 2
dim = [];
end
end
end
% Check that either dim or box is provided
if (isempty (dim)) && (isempty (box))
stk_error (['The dimension argument can be omitted if, and only if, a ' ...
'valid box argument is provided instead.'], 'InvalidArgument');
end
% Process box argument
if isempty (box)
colnames = {};
else
box = stk_hrect (box); % convert input argument to a proper box
colnames = box.colnames;
if isempty (dim)
dim = size (box, 2);
elseif dim ~= size (box, 2)
stk_error (['The dimension argument must be compatible with' ...
'the box argument when both are provided.'], 'InvalidArgument');
end
end
% Generate a Sobol sequence in [0; 1]^dim
x = transpose (__stk_sampling_sobol_mex__ (n, dim, do_skip));
% Create dataframe output
x = stk_dataframe (x, colnames);
if ~ isempty (box),
x = stk_rescale (x, [], box);
end
end % function
|
