# Question:How to make this procedure more efficient?

## Question:How to make this procedure more efficient?

Maple

Hi,

This thread is more or less related to a previous one about the Statistics:-Sample procedure.
(see https://www.mapleprimes.com/questions/228421-A-Serious-Problem-With-StatisticsSample )
I've just implemented two variant of the Box-Muller procedure to sample normal rvs.
The source is "The art of computer programming", Donald E. Knuth, 2nd edition, p117 (aka "algorithm P").

The first implementation (BoxMuller_1) is basically what D.E. Knuth writes, except that I "vectorize" some operations in order to avoid using an if.. then..else structure (as a minor consequence I generate a little bit more numbers than required).
This procedure uses the build-in Maple's procedure select (please see the link above and acer's and carl love's replies).
It appears to be relatively slow.
More of this CodeTools:-Usage(BoxMuller_1(10^6)) generates a "conenction to kernel lost" for some unknown reason.

The variant named BoxMuller_2 uses sort and ListTools:-BinaryPlace instead of select.
Ir appears to be around 3 times faster than BoxMuller_1, but remains 10 time slower than Statistics:-Sample(..., method=envelope) (here again, see the link above to understand why method=envelope is needed).

I wonder how it could be possible to speed up this procedure. In particular acer showed in one of his reply (again see the link above) how using hfloats can improve the efficiency of a procedure, but I'm very incompetent on this point.
Does anyone have any idea?

PS: this file has been written with Maple 2015.2

 > restart:
 > BoxMuller_1 := proc(N)   local V, S, T, L, X1, X2:   V  := Statistics:-Sample(Uniform(-1, 1), [ceil(2*N/3), 2]):   S  := V[..,1]^~2 +~ V[..,2]^~2;   T  := < subs(NULL=infinity, select(`<`, S, 1)) | V >;   T[.., 1]  := (-2 *~ log~(T[.., 1]) /~ T[.., 1])^~(0.5);   X1 := select[flatten](type, T[.., 1] *~  T[..,2], 'float');   X2 := select[flatten](type, T[.., 1] *~  T[..,3], 'float');   return : end proc: BoxMuller_2 := proc(N)   local V1, V2, S, W, u, r, T, X1, X2:   V1 := Statistics:-Sample(Uniform(-1, 1), ceil(2*N/3)):   V2 := Statistics:-Sample(Uniform(-1, 1), ceil(2*N/3)):   S  := V1^~2 +~ V2^~2;   W  := sort(S, output = [sorted, permutation]):   u  := ListTools:-BinaryPlace(W[1], 1);   r  := [\$1..u]:   T  := (-2 *~ log~(W[1][r]) /~ W[1][r])^~(0.5);   X1 := T *~ V1[W[2][r]];   X2 := T *~ V2[W[2][r]];   return ^+: end proc:
 > CodeTools:-Usage(Statistics:-Sample(Uniform(-1, 1), 10^5)): CodeTools:-Usage(Statistics:-Sample(Uniform(-1, 1), 10^5, method=envelope)): # CodeTools:-Usage(BoxMuller_1(10^6)): # generates a "connection to kernel lost" error msg CodeTools:-Usage(BoxMuller_2(10^6)): print(): CodeTools:-Usage(BoxMuller_1(10^5)): CodeTools:-Usage(BoxMuller_2(10^5)): print(): S1 := BoxMuller_1(10^5): S2 := BoxMuller_2(10^5):
 memory used=1.17MiB, alloc change=0 bytes, cpu time=12.00ms, real time=12.00ms, gc time=0ns memory used=3.32MiB, alloc change=32.00MiB, cpu time=64.00ms, real time=64.00ms, gc time=0ns memory used=148.17MiB, alloc change=119.93MiB, cpu time=706.00ms, real time=637.00ms, gc time=14.63ms
 memory used=67.43MiB, alloc change=186.80MiB, cpu time=946.00ms, real time=762.00ms, gc time=285.14ms memory used=14.81MiB, alloc change=0 bytes, cpu time=67.00ms, real time=62.00ms, gc time=0ns
 (1)
 > if false then DocumentTools:-Tabulate(   [     plots:-display(       Statistics:-Histogram(S1),       plot(Statistics:-PDF(Normal(0, 1), x), x=-4..4, color=red,thickness=3)     ),     plots:-display(       Statistics:-Histogram(S2),       plot(Statistics:-PDF(Normal(0, 1), x), x=-4..4, color=red,thickness=3)     )   ], width=60 ) end if:
 > # carl love's procedure with slight modifications t SampleCheck := proc(X, f, N::posint) uses St= Statistics;   local S, M, O, E:   S:= f(N):   M:= numelems(S):   O:= )>:   E:= St:-Probability~(X <~ (rhs@lhs)~(O), 'numeric') * M:    end proc:
 > interface(rtablesize=20): SampleCheck(Statistics:-RandomVariable(Normal(0,1)), 'BoxMuller_1', 10^5), SampleCheck(Statistics:-RandomVariable(Normal(0,1)), 'BoxMuller_2', 10^5);
 (2)
 >