qmc_sequence

@info:
! @deftypefn {Mex File} {[@var{a}, @var{s}, @var{info}] =} qmc_sequence (@var{d}, @var{s}, @var{t}, @var{n}, @var{lu})
! @anchor{qmc_sequence}
! @sp 1
! Computes quasi Monte-Carlo sequence (Sobol numbers).
! @sp 2
! @strong{Inputs}
! @sp 1
! @table @ @var
! @item d
! Integer scalar, dimension.
! @item s
! Integer scalar (int64), seed.
! @item t
! Integer scalar, sequence type:
!  @sp 1
!  @table @ @samp
!  @item @var{t}=0
!  Uniform numbers in a n-dimensional (unit by default) hypercube 
!  @item @var{t}=1
!  Gaussian numbers
!  @item @var{t}=2
!  Uniform numbers on a n-dimensional (unit by default) hypersphere
!  @end table
! @item n
! Integer scalar, number of elements in the sequence.
! @item lu
! Optional argument, the interpretation depends on its size:
!  @sp 1
!  @table @ @samp
!  @item @var{d}x2 array of doubles
!  Lower and upper bounds of the hypercube (default is 0-1 in all dimensions). @var{t} must be equal to zero.
!  @item @var{d}x@var{d} array of doubles
!  Lower cholesky of the covariance matrix of the Gaussian variates (default is the identity matrix). @var{t} must be equal to one.
!  @item scalar double
!  Radius of the hypershere (default is one). @var{t} must be equal to two.
!  @end table
! @end table
! @sp 2
! @strong{Outputs}
! @sp 1
! @table @ @var
! @item a
! @var{n}x@var{d} matrix of doubles, the Sobol sequence.
! @item s
! Integer scalar (int64), current value of the seed.
! @item info
! Integer scalar, equal to 1 if mex routine fails, 0 otherwise.
! @end table
! @sp 2
! @strong{This function is called by:}
! @sp 2
! @strong{This function calls:}
! @sp 2
! @end deftypefn
@eod:

0001 %@info:
0002 %! @deftypefn {Mex File} {[@var{a}, @var{s}, @var{info}] =} qmc_sequence (@var{d}, @var{s}, @var{t}, @var{n}, @var{lu})
0003 %! @anchor{qmc_sequence}
0004 %! @sp 1
0005 %! Computes quasi Monte-Carlo sequence (Sobol numbers).
0006 %! @sp 2
0007 %! @strong{Inputs}
0008 %! @sp 1
0009 %! @table @ @var
0010 %! @item d
0011 %! Integer scalar, dimension.
0012 %! @item s
0013 %! Integer scalar (int64), seed.
0014 %! @item t
0015 %! Integer scalar, sequence type:
0016 %!  @sp 1
0017 %!  @table @ @samp
0018 %!  @item @var{t}=0
0019 %!  Uniform numbers in a n-dimensional (unit by default) hypercube
0020 %!  @item @var{t}=1
0021 %!  Gaussian numbers
0022 %!  @item @var{t}=2
0023 %!  Uniform numbers on a n-dimensional (unit by default) hypersphere
0024 %!  @end table
0025 %! @item n
0026 %! Integer scalar, number of elements in the sequence.
0027 %! @item lu
0028 %! Optional argument, the interpretation depends on its size:
0029 %!  @sp 1
0030 %!  @table @ @samp
0031 %!  @item @var{d}x2 array of doubles
0032 %!  Lower and upper bounds of the hypercube (default is 0-1 in all dimensions). @var{t} must be equal to zero.
0033 %!  @item @var{d}x@var{d} array of doubles
0034 %!  Lower cholesky of the covariance matrix of the Gaussian variates (default is the identity matrix). @var{t} must be equal to one.
0035 %!  @item scalar double
0036 %!  Radius of the hypershere (default is one). @var{t} must be equal to two.
0037 %!  @end table
0038 %! @end table
0039 %! @sp 2
0040 %! @strong{Outputs}
0041 %! @sp 1
0042 %! @table @ @var
0043 %! @item a
0044 %! @var{n}x@var{d} matrix of doubles, the Sobol sequence.
0045 %! @item s
0046 %! Integer scalar (int64), current value of the seed.
0047 %! @item info
0048 %! Integer scalar, equal to 1 if mex routine fails, 0 otherwise.
0049 %! @end table
0050 %! @sp 2
0051 %! @strong{This function is called by:}
0052 %! @sp 2
0053 %! @strong{This function calls:}
0054 %! @sp 2
0055 %! @end deftypefn
0056 %@eod:
0057 
0058 %@test:1
0059 %$ t = ones(3,1);
0060 %$
0061 %$ d = 2;
0062 %$ n = 100;
0063 %$ s = int64(0);
0064 %$
0065 %$ try
0066 %$   [draws, S] = qmc_sequence(d,s,0,n);
0067 %$ catch
0068 %$   t(1) = 0;
0069 %$ end
0070 %$
0071 %$ try
0072 %$   [draws, S] = qmc_sequence(d,s,1,n);
0073 %$ catch
0074 %$   t(2) = 0;
0075 %$ end
0076 %$
0077 %$ try
0078 %$   [draws, S] = qmc_sequence(d,s,2,n);
0079 %$ catch
0080 %$   t(3) = 0;
0081 %$ end
0082 %$
0083 %$ T = all(t);
0084 %@eof:1
0085 
0086 %@test:2
0087 %$ t = ones(3,1);
0088 %$
0089 %$ d = 2;
0090 %$ n = 100;
0091 %$ s = int64(0);
0092 %$
0093 %$ [draws1, S] = qmc_sequence(d,s,0,n);
0094 %$ [draws2, Q] = qmc_sequence(d,S,0,n);
0095 %$ [draws3, P] = qmc_sequence(d,s,0,2*n);
0096 %$
0097 %$ t(1) = dyn_assert(s,int64(0));
0098 %$ t(2) = dyn_assert(P,Q);
0099 %$ t(3) = dyn_assert([draws1,draws2],draws3);
0100 %$ T = all(t);
0101 %@eof:2
0102 
0103 %@test:3
0104 %$
0105 %$ d = 2;
0106 %$ n = 100;
0107 %$ s = int64(0);
0108 %$
0109 %$ [draws1, S] = qmc_sequence(d,s,0,n,[0 , 2; -1, 2]);
0110 %$ [draws2, Q] = qmc_sequence(d,s,0,n);
0111 %$
0112 %$ draws3 = draws2;
0113 %$ draws3(1,:) = 2*draws2(1,:);
0114 %$ draws3(2,:) = 3*draws2(2,:)-1;
0115 %$ t(1) = dyn_assert(S,Q);
0116 %$ t(2) = dyn_assert(draws1,draws3);
0117 %$ T = all(t);
0118 %@eof:3
0119 
0120 %@test:4
0121 %$
0122 %$ d = 2;
0123 %$ n = 100;
0124 %$ s = int64(0);
0125 %$ radius = pi;
0126 %$
0127 %$ [draws, S] = qmc_sequence(d,s,2,n,radius);
0128 %$
0129 %$ t(1) = dyn_assert(sqrt(draws(:,3)'*draws(:,3)),radius,1e-14);
0130 %$ t(2) = dyn_assert(sqrt(draws(:,5)'*draws(:,5)),radius,1e-14);
0131 %$ t(3) = dyn_assert(sqrt(draws(:,7)'*draws(:,7)),radius,1e-14);
0132 %$ t(4) = dyn_assert(sqrt(draws(:,11)'*draws(:,11)),radius,1e-14);
0133 %$ t(5) = dyn_assert(sqrt(draws(:,13)'*draws(:,13)),radius,1e-14);
0134 %$ t(6) = dyn_assert(sqrt(draws(:,17)'*draws(:,17)),radius,1e-14);
0135 %$ t(7) = dyn_assert(sqrt(draws(:,19)'*draws(:,19)),radius,1e-14);
0136 %$ t(8) = dyn_assert(sqrt(draws(:,23)'*draws(:,23)),radius,1e-14);
0137 %$ t(9) = dyn_assert(sqrt(draws(:,29)'*draws(:,29)),radius,1e-14);
0138 %$ T = all(t);
0139 %@eof:4
0140 
0141 %@test:5
0142 %$
0143 %$ d = 2;
0144 %$ n = 100000;
0145 %$ b = 100;
0146 %$ s = int64(5);
0147 %$
0148 %$ covariance = [.4 -.1; -.1 .2];
0149 %$ chol_covariance = transpose(chol(covariance));
0150 %$
0151 %$ draws = [];
0152 %$
0153 %$ for i=1:b
0154 %$     [tmp, s] = qmc_sequence(d,s,1,n,chol_covariance);
0155 %$     draws = [draws, tmp];
0156 %$ end
0157 %$
0158 %$ COVARIANCE = draws*draws'/(b*n);
0159 %$
0160 %$ t(1) = dyn_assert(covariance,COVARIANCE,1e-6);
0161 %$ T = all(t);
0162 %@eof:5