# Question:How to fasten this simulation?

## Question:How to fasten this simulation?

Maple 2015

Hi,

I want to simulate a decision rule (let's say a production control strategy to fix the ideas) based on the outcomes of a binomial random variable.

One stage of this simulation involves operations which correspond to the NaiveProc described below. This procedure is a straightforward and rather naive translation in Maple of these operations.
It is given here to help you to understand what is to be done.

This coding being extremely slow, I implemented the required operations differently.
Given the fact that I'm only concerned by binomial samplings with outcomes 0 and 1, I basically replaced these samplings by three operations  U:=... , Got01:=... , Got0:= ... .
Each new implementation is about 200 times faster than the naive one.

My question is: Can we do still better?
For information the same simulation realized with R runs in less than 15s when REP = 10^8 (more than 100 times faster than what I'm capable to do with Maple)

 > restart:
 > with(Statistics):
 > # This is the naive coding of the problem # # It has the advantage to explain clearly (at least I hope so) what is to be done. NaiveProc := proc(REP, a, b, N)   local roll, rep, p, K, N0, N1:   roll := rand(a .. b);   N0   := 0:   N1   := 0:   for rep from 1 to REP do     p := roll();     K := Sample(Binomial(N, p), 1)[1]:     if K=0 then       N0 := N0+1:     elif K=1 then       N1 := N1+1:     end if:   end do:   N0, N1: end proc:      CodeTools:-Usage(NaiveProc(10^3, 0, 5e-4, 25000));
 memory used=78.92MiB, alloc change=68.00MiB, cpu time=1.21s, real time=6.30s, gc time=24.80ms
 (1)
 > # As the previous code is very slow I implemented several variants based on the fact # that I care only on values of K equal to 0 and 1 and that it is then useless to sample # a Binomial random variable to do this. # # Each of them is much faster than the naive coding... but can we even faster? # For comparison, the same simulation realized in R takes less than 15s to run with REP=10^8 pi := (n, p, k) -> binomial(n, k) * (1-p)^(n-k) * p^k; f1 := proc(REP, a, b, N)   local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:   Q     := Sample(Uniform(a, b), REP)^+:   P0    := pi~(N, Q, 0):   P1    := pi~(N, Q, 1):   P01   := P0 + P1:   U     := Sample(Uniform(0, 1), REP)^+:   Got01 := Select(t -> is(t <= 0), U-P01):   Got0  := Select(t -> is(t <= 0), U-P0 );   N01   := numelems(Got01):   N0    := numelems(Got0):   N1    := N01-N0:   N0, N1; end proc: CodeTools:-Usage(f1(10^5, 0, 5e-4, 25000)); print(): f2 := proc(REP, a, b, N)   local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:   Q     := Sample(Uniform(a, b), REP)^+:   P0    := pi~(N, Q, 0):   P1    := pi~(N, Q, 1):   P01   := P0 + P1:   U     := Sample(Uniform(0, 1), REP)^+:   Got01 := select[flatten](`<=`, U-P01, 0):   Got0  := select[flatten](`<=`, U-P0 , 0):   N01   := numelems(Got01):   N0    := numelems(Got0):   N1    := N01-N0:   N0, N1; end proc: CodeTools:-Usage(f2(10^5, 0, 5e-4, 25000)); print(): f3 := proc(REP, a, b, N)   local Q, P0, P1, P01, U, Got01, Got0, N01, N0, N1:   Q     := Sample(Uniform(a, b), REP)^+:   P0    := pi~(N, Q, 0):   P1    := pi~(N, Q, 1):   P01   := P0 + P1:   U     := Sample(Uniform(0, 1), REP)^+:   Got01 := Vector(REP, i -> `if`(U[i] <= P01[i], 1, 0)):   Got0  := Vector(REP, i -> `if`(U[i] <= P0 [i], 1, 0)):   N01   := add(Got01):   N0    := add(Got0):   N1    := N01-N0:   N0, N1; end proc: CodeTools:-Usage(f3(10^5, 0, 5e-4, 25000));
 memory used=376.83MiB, alloc change=296.77MiB, cpu time=4.12s, real time=4.03s, gc time=430.38ms
 memory used=173.26MiB, alloc change=0 bytes, cpu time=2.30s, real time=1.98s, gc time=461.22ms
 memory used=154.19MiB, alloc change=0 bytes, cpu time=2.28s, real time=2.01s, gc time=432.40ms
 (2)
 >
 >