It depends on the distribution. For some distributions (e.g. normal) a specialized algorithm is used. For normal this is the ziggurat method. In the generic case, the method is adaptive rejection applied to the log transformed density or in some cases to the original density. It is decided by some heuristics embedded in the code.
One possibility is to perform all of the the preprocessing once. This should make things a bit faster. Consider:
restart:
Digits :=10:
test3:=proc(N,n)
local i,a,X;
use Statistics in
X := RandomVariable(Normal(0,1)):
for i from 1 to N do
a:=Sample(X,n):
end do:
end use;
end proc:
time(test3(30000,100));
12.733
Note, that the call to Sample below creates a procedure for generating random samples drawn from the normal distribution. So it is not longer necessary to do the parameter processing over and over again.
restart:
Digits :=10:
test3:=proc(N,n)
local i,a,X, S;
use Statistics in
X := RandomVariable(Normal(0,1)): S := Sample(X);
for i from 1 to N do
a:=S(n):
end do:
end use;
end proc:
time(test3(30000,100));
1.446
Of course, it is even more efficient to do the simulation once.
restart:
Digits :=10:
test3:=proc(N,n)
local i,a,X, S;
use Statistics in
X := RandomVariable(Normal(0,1)): S := Sample(X);
a := ArrayTools:-Alias(S(N*n), [1..N, 1..n]);
end use;
end proc:
time(test3(30000,100));
0.268
restart:
Digits :=10;
test4:=proc(N,n)
local i,a,b,X;
use Statistics in
X := RandomVariable(Normal(0,1)):
a:=Sample(X,n):
for i from 1 to N do
b:=Mean(a)
end do:
end use;
end proc:
time(test4(30000,100));
10
7.126
Again, we can create a procedure for computing the mean of a data sample.
restart:
Digits :=10;
test4:=proc(N,n)
local i,a,b,X,M;
use Statistics in
X := RandomVariable(Normal(0,1)):
a:=Sample(X,n): M := MakeProcedure(Mean, sample);
for i from 1 to N do
b:=M(a);
end do:
end use;
end proc:
time(test4(30000,100));
10
0.277
Try running the command-line version (maple instead of xmaple).
-Alex
The following method (due to Abate and Valko) works pretty well in most cases. One of it's most interesting features is that going to high precision allows you to speed up computations by quite a bit. Note that this is a very rough prototype.
Talbot := proc(f, t, M)
local r, ss, ds, F0, F1, F;
Digits := M;
r := .4*M/t;
ss := theta -> r*theta*(cot(theta)+I):
ds := theta -> theta+(theta*cot(theta)-1)*cot(theta):
F0 := .5*f(r)*exp(r*t);
F1 := add(evalf('Re'(exp(t*ss(k*Pi/M))*f(ss(k*Pi/M))*(1+I*ds(k*Pi/M)))), k = 1..M-1):
return evalf(r/M*(F0+F1));
end proc: # Talbot
> f := t -> sin(sqrt(t));
f := t -> sin(sqrt(t))
> ff := unapply(inttrans[laplace](f(t), t, u), u);
1/2
Pi exp(-1/4 1/u)
ff := u -> 1/2 -------------------
3/2
u
> evalf(f(.17) = Talbot(ff, .17, 20));
0.4007273271 = 0.4007273271
> evalf(f(7) = Talbot(ff, 7, 20));
0.4757718382 = 0.4757718382
-Alex Potapchik
Correct me if I am wrong, but I think `i` has to be declared local to CS. This will change the timings. :(
-Alex
BM := proc(n::posint, m::posint)
local X;
uses Statistics,ListTools;
X := RandomVariable(Normal(0,1/sqrt(n)));
'<0,PartialSums([seq(i,i=Sample(X,n))])[]>'$m;
end proc:
BM2 := proc(n::posint, m::posint)
local A, i;
A := ArrayTools:-Alias(Statistics:-Sample(Normal(0, 1/sqrt(n)), m*n), [1..m, 1..n]);
for i to m do
A[i] := Statistics:-CumulativeSum(A[i]);
end do;
A;
end proc;
> time(BM(200, 200));
0.924
> time(BM2(200, 200));
0.048
> time(BM2(1000, 1000));
0.956
-Alex
> restart; f := simplify(sqrt(1/x^4+1)) assuming 1 <= x, x <= 2;
4 1/2
(1 + x )
f := -----------
2
x
> int(f, x = 1..2);
1/2 1/2 1/2 1/2 1/2
3 17 2 2 2 2
------- + EllipticK(----) - EllipticF(4/5, ----) - 2 EllipticE(----) + 2 EllipticE(4/5, ----)
10 2 2 2 2
> evalf(%);
1.132090394
> evalf(Int(f, x = 1..2));
1.132090393
- Alex
Here is a simple procedure which would sample from a multinormal distribution with the specified covariance matrix:
> MultivariateNormalSample := proc(Sigma, N)
> local A, R, d;
> d := LinearAlgebra[RowDimension](Sigma);
> R := Matrix(LinearAlgebra[LUDecomposition](evalf(Sigma), 'method' = 'Cholesky'), datatype = float[8]);
> A := Statistics[Sample](Normal(0, 1), d*N);
> A := ArrayTools[Alias](A, [1..d, 1..N]);
> rtable_options(A, subtype = Matrix);
> return R.A;
> end proc:
Let's try.
> Sigma := <<1., .3>|<.3, 1.>>;
[1. 0.3]
Sigma := [ ]
[0.3 1. ]
> S := MultivariateNormalSample(Sigma, 10^5);
[ 2 x 100000 Matrix ]
S := [ Data Type: float[8] ]
[ Storage: rectangular ]
[ Order: Fortran_order ]
> map(evalf[5], Statistics[CovarianceMatrix](LinearAlgebra[Transpose](S)));
[0.99555 0.29937]
[ ]
[0.29937 0.99478]
> Sigma := <<2.0, 0.3, 0.2>|<0.3, 1.0, 0.5>|<0.2, 0.5, 3.0>>;
[2.0 0.3 0.2]
[ ]
Sigma := [0.3 1.0 0.5]
[ ]
[0.2 0.5 3.0]
> S := MultivariateNormalSample(Sigma, 10^5);
[ 3 x 100000 Matrix ]
S := [ Data Type: float[8] ]
[ Storage: rectangular ]
[ Order: Fortran_order ]
> map(evalf[5], Statistics[CovarianceMatrix](LinearAlgebra[Transpose](S)));
[2.0111 0.30239 0.20690]
[ ]
[0.30239 0.99475 0.50142]
[ ]
[0.20690 0.50142 3.0098 ]
-Alex, Maplesoft
Here is a simple procedure which would sample from a multinormal distribution with the specified covariance matrix:
> MultivariateNormalSample := proc(Sigma, N)
> local A, R, d;
> d := LinearAlgebra[RowDimension](Sigma);
> R := Matrix(LinearAlgebra[LUDecomposition](evalf(Sigma), 'method' = 'Cholesky'), datatype = float[8]);
> A := Statistics[Sample](Normal(0, 1), d*N);
> A := ArrayTools[Alias](A, [1..d, 1..N]);
> rtable_options(A, subtype = Matrix);
> return R.A;
> end proc:
Let's try.
> Sigma := <<1., .3>|<.3, 1.>>;
[1. 0.3]
Sigma := [ ]
[0.3 1. ]
> S := MultivariateNormalSample(Sigma, 10^5);
[ 2 x 100000 Matrix ]
S := [ Data Type: float[8] ]
[ Storage: rectangular ]
[ Order: Fortran_order ]
> map(evalf[5], Statistics[CovarianceMatrix](LinearAlgebra[Transpose](S)));
[0.99555 0.29937]
[ ]
[0.29937 0.99478]
> Sigma := <<2.0, 0.3, 0.2>|<0.3, 1.0, 0.5>|<0.2, 0.5, 3.0>>;
[2.0 0.3 0.2]
[ ]
Sigma := [0.3 1.0 0.5]
[ ]
[0.2 0.5 3.0]
> S := MultivariateNormalSample(Sigma, 10^5);
[ 3 x 100000 Matrix ]
S := [ Data Type: float[8] ]
[ Storage: rectangular ]
[ Order: Fortran_order ]
> map(evalf[5], Statistics[CovarianceMatrix](LinearAlgebra[Transpose](S)));
[2.0111 0.30239 0.20690]
[ ]
[0.30239 0.99475 0.50142]
[ ]
[0.20690 0.50142 3.0098 ]
-Alex, Maplesoft
You should be able to use this library:
http://parallel.bas.bg/~emanouil/sequences.html.
This is probably the fastest implementation of sobol sequences available at this moment. A lot of it is assembli code. Unfortunately it only works with gcc.
-Alex, Maplesoft
Right. The function has to be exported. Sorry I forgot to mention this. I used the same method as suggested by Alec. The watcom compiler also provides options for specifying which functions have to be exported.
-Alex
Right. The function has to be exported. Sorry I forgot to mention this. I used the same method as suggested by Alec. The watcom compiler also provides options for specifying which functions have to be exported.
-Alex
Note that this can break other things in your maple session. Consider:
> restart; `property/ConvertRelation`:=()->FAIL:
> is(x > 1) assuming x > 2;
Error, (in assuming) when calling 't1'. Received: 'attempting to assign to `FAIL` which is protected'
Can you use unevaluated piecewise or another data structure, which prints the same way?
> vf := 3*vx[1] < [6, 1] and vx[1] < [-2\
> , 5], [3, 4]+3*vx[1], 2*vx[1] < [8, -4\
> ] and -vx[1] < [2, -5], [1, 9]+2*vx[1]\
> , -3*vx[1] < [-6, -1] and -2*vx[1] < [-8, 4], [9, 5];
vf := 3 vx[1] < [6, 1] and vx[1] < [-2, 5], [3, 4] + 3 vx[1],
2 vx[1] < [8, -4] and -vx[1] < [2, -5], [1, 9] + 2 vx[1],
-3 vx[1] < [-6, -1] and -2 vx[1] < [-8, 4], [9, 5]
> 'piecewise'(vf);
{ [3, 4] + 3 vx[1] 3 vx[1] < [6, 1] and vx[1] < [-2, 5]
{
{ [1, 9] + 2 vx[1] 2 vx[1] < [8, -4] and -vx[1] < [2, -5]
{
{ [9, 5] -3 vx[1] < [-6, -1] and -2 vx[1] < [-8, 4]
> %; # This will not work
Error, (in property/ConvertRelation) invalid relation arguments
> `print/mypiecewise` := eval(`print/piecewise`);
print/mypiecewise := proc() ... end proc
> mypiecewise(vf);
{ [3, 4] + 3 vx[1] 3 vx[1] < [6, 1] and vx[1] < [-2, 5]
{
{ [1, 9] + 2 vx[1] 2 vx[1] < [8, -4] and -vx[1] < [2, -5]
{
{ [9, 5] -3 vx[1] < [-6, -1] and -2 vx[1] < [-8, 4]
Note that this can break other things in your maple session. Consider:
> restart; `property/ConvertRelation`:=()->FAIL:
> is(x > 1) assuming x > 2;
Error, (in assuming) when calling 't1'. Received: 'attempting to assign to `FAIL` which is protected'
Can you use unevaluated piecewise or another data structure, which prints the same way?
> vf := 3*vx[1] < [6, 1] and vx[1] < [-2\
> , 5], [3, 4]+3*vx[1], 2*vx[1] < [8, -4\
> ] and -vx[1] < [2, -5], [1, 9]+2*vx[1]\
> , -3*vx[1] < [-6, -1] and -2*vx[1] < [-8, 4], [9, 5];
vf := 3 vx[1] < [6, 1] and vx[1] < [-2, 5], [3, 4] + 3 vx[1],
2 vx[1] < [8, -4] and -vx[1] < [2, -5], [1, 9] + 2 vx[1],
-3 vx[1] < [-6, -1] and -2 vx[1] < [-8, 4], [9, 5]
> 'piecewise'(vf);
{ [3, 4] + 3 vx[1] 3 vx[1] < [6, 1] and vx[1] < [-2, 5]
{
{ [1, 9] + 2 vx[1] 2 vx[1] < [8, -4] and -vx[1] < [2, -5]
{
{ [9, 5] -3 vx[1] < [-6, -1] and -2 vx[1] < [-8, 4]
> %; # This will not work
Error, (in property/ConvertRelation) invalid relation arguments
> `print/mypiecewise` := eval(`print/piecewise`);
print/mypiecewise := proc() ... end proc
> mypiecewise(vf);
{ [3, 4] + 3 vx[1] 3 vx[1] < [6, 1] and vx[1] < [-2, 5]
{
{ [1, 9] + 2 vx[1] 2 vx[1] < [8, -4] and -vx[1] < [2, -5]
{
{ [9, 5] -3 vx[1] < [-6, -1] and -2 vx[1] < [-8, 4]
Another option to consider is calling your compiled C routine directly from Maple (see ?define_external). This way you will be able to run your simulation as efficiently as possible and at the same time you will be able to use the results of your simulation in a Maple session to do any further analysis.
