## 16 Reputation

15 years, 154 days

## It depends on the...

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...

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

## This may be UI-related...

Try running the command-line version (maple instead of xmaple). -Alex

## The following method (due to...

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,...

Correct me if I am wrong, but I think `i` has to be declared local to CS. This will change the timings. :( -Alex

## Here is a slightly more efficient way to...

```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

## Possible workaround....

```> 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...

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);
>     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  ]
[ Storage: rectangular ]
[ Order: Fortran_order ]

> map(evalf, 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  ]
[ Storage: rectangular ]
[ Order: Fortran_order ]

> map(evalf, 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...

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);
>     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  ]
[ Storage: rectangular ]
[ Order: Fortran_order ]

> map(evalf, 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  ]
[ Storage: rectangular ]
[ Order: Fortran_order ]

> map(evalf, 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 can use this...

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...

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...

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...

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 < [6, 1] and vx < [-2\ > , 5], [3, 4]+3*vx, 2*vx < [8, -4\ > ] and -vx < [2, -5], [1, 9]+2*vx\ > , -3*vx < [-6, -1] and -2*vx < [-8, 4], [9, 5]; vf := 3 vx < [6, 1] and vx < [-2, 5], [3, 4] + 3 vx, 2 vx < [8, -4] and -vx < [2, -5], [1, 9] + 2 vx, -3 vx < [-6, -1] and -2 vx < [-8, 4], [9, 5] > 'piecewise'(vf); { [3, 4] + 3 vx 3 vx < [6, 1] and vx < [-2, 5] { { [1, 9] + 2 vx 2 vx < [8, -4] and -vx < [2, -5] { { [9, 5] -3 vx < [-6, -1] and -2 vx < [-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 3 vx < [6, 1] and vx < [-2, 5] { { [1, 9] + 2 vx 2 vx < [8, -4] and -vx < [2, -5] { { [9, 5] -3 vx < [-6, -1] and -2 vx < [-8, 4] ```

## Note that this can break...

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 < [6, 1] and vx < [-2\ > , 5], [3, 4]+3*vx, 2*vx < [8, -4\ > ] and -vx < [2, -5], [1, 9]+2*vx\ > , -3*vx < [-6, -1] and -2*vx < [-8, 4], [9, 5]; vf := 3 vx < [6, 1] and vx < [-2, 5], [3, 4] + 3 vx, 2 vx < [8, -4] and -vx < [2, -5], [1, 9] + 2 vx, -3 vx < [-6, -1] and -2 vx < [-8, 4], [9, 5] > 'piecewise'(vf); { [3, 4] + 3 vx 3 vx < [6, 1] and vx < [-2, 5] { { [1, 9] + 2 vx 2 vx < [8, -4] and -vx < [2, -5] { { [9, 5] -3 vx < [-6, -1] and -2 vx < [-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 3 vx < [6, 1] and vx < [-2, 5] { { [1, 9] + 2 vx 2 vx < [8, -4] and -vx < [2, -5] { { [9, 5] -3 vx < [-6, -1] and -2 vx < [-8, 4] ```

## External calling...

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.
 1 2 Page 1 of 2
﻿