# How can I improve performance of my simulation code?

1) Generate nNsamples from beta distribution (a, b) and the parameter of wighted c

Subscript[a, 0] = 2;
Subscript[b, 0] = 4;
nN = 1000;
n = 50;
α = 0.05;
n1 = nN*α;
c = 2;

w1 =
Table[
{i, nN}, {j, 1, n}];

me = Subscript[a, 0]/(Subscript[a, 0] + Subscript[b, 0])

se =
{(Subscript[a, 0] Subscript[b, 0])/
((Subscript[a, 0] + Subscript[b, 0])^2 (Subscript[a, 0] + Subscript[b, 0] + 1))}


2) Calculate unknown parameters a,b for each random sample*)

m6 = Table[Mean[w1[[i, All]]], {i, nN}];

s6 = Table[Variance[w1[[i, All]]], {i, nN}];
e6 = Table[m6[[i]] - a/(a + b) == 0, {i, nN}];
e26 = Table[s6[[i]] - (a b)/((a + b)^2 (a + b + 1)) == 0, {i, nN}];

v1 =
Table[
FindRoot[{e6[[i]], e26[[i]]}, {{a, 1.5}, {b, 3.5}}, MaxIterations -> 5],
{i, nN}];}

OverHat[a] = a /. v1;
OverHat[b] = b /. v1;
f1 = Mean[OverHat[a]]
f2 = Mean[OverHat[b]]


3) Calculate suggested test statistical

T =
Table[
(1/n)*
Sum[
(w1[[i, l]]*w1[[i, k]] - 2*(w1[[i, l]]*m6[[i]]) + m6[[i]]^2)/
((w1[[i, l]] + w1[[i, k]] + c)*m6[[i]]^2) +
(2*(w1[[i, l]]*w1[[i, k]]*(1 - w1[[i, k]]) -
(w1[[i, l]]*m6[[i]])*(1 - w1[[i, l]])))/
((w1[[i, l]] + w1[[i, k]] + c)^2*OverHat[a][[i]]*m6[[i]]) +
(2*(w1[[i, l]]*w1[[i, k]]*(w1[[i, k]]*(w1[[i, l]] - 2) + 1)))/
((w1[[i, l]] + w1[[i, k]] + c)^3*OverHat[a][[i]]^2),
{l, 1, n}, {k, 1, n}],
{i, nN}];


4) Simulate (n1) bootstrap samples from a beta distribution with parameters OverHat[a], OverHat[b])

w2 =
Table[
{i, nN}, {j, 1, n}];
w3 = Table[RandomChoice[w2[[i]], {n1, n}], {i, nN}];

Dimensions[w3]

m7 = Table[Mean[w3[[i, r, All]]], {i, nN}, {r, n1}];

s7 = Table[Variance[w3[[i, r, All]]], {i, nN}, {r, n1}];}


5) Calculate for each bootstrap sample the value of the test statistic

v2 =
Table[
FindRoot[
{m7[[i, r]] - a1/(a1 + b1) == 0,
s7[[i, r]] - (a1 b1)/((a1 + b1)^2 (a1 + b1 + 1)) == 0},
{a1, f1}, {b1, f2},
MaxIterations -> 5],
{i, nN}, {r, n1}];

OverHat[a1] = a1 /. v2;

SuperStar[T] =
(Table[(1/n)*
Sum[
(w3[[i, k, l]]*w3[[i, k, z]] - 2*(w3[[i, k, l]]*m7[[i, k]]) + m7[[i, k]]^2)/
((w3[[i, k, l]] + w3[[i, k, z]] + c)*m7[[i, k]]^2) +
(2*(w3[[i, k, l]]*w3[[i, k, z]]*(1 - w3[[i, k, z]]) -
(w3[[i, k, l]]*m7[[i, k]])*(1 - w3[[i, k, l]])))/
((w3[[i, k, l]] + w3[[i, k, z]] + c)^2*OverHat[a1][[i, k]]*m7[[i, k]]) +
(2*(w3[[i, k, l]]*w3[[i, k, z]]*(w3[[i, k, z]]*(w3[[i, k, l]] - 2) + 1)))/
((w3[[i, k, l]] + w3[[i, k, z]] + c)^3*OverHat[a1][[i, k]]^2),
{l, 1, n}, {z, 1, n}],
{i, nN}, {k, n1}];)

Dimensions[SuperStar[T]]}


6) Calculate p_value as order statistic of the bootstrap sample

arr = Sort /@ SuperStar[T];
Dimensions[arr]
d = n1 - (n1*α)
l = Round[d]
pn = arr[[All, l]];
Dimensions[pn]}


7) Calculate power of test (Reject the hypothesis if the test statistic

T of your original sample is greater than p_n)

pp = Total[Positive[T - pn]]
repp = Total[UnitStep[T - pn]]
repp2 = Total[UnitStep[pn - T]]}

• so whats your question? First you need to find out which step takes the most time and then you can maybe ask a more specific question. Try AbsoluteTime[] and have a look at mathematica.stackexchange.com/questions/29349/… Commented Oct 15, 2017 at 13:16
• in my simulation step 5 takes the most time,precisely when calculate SuperStar[T] Commented Oct 15, 2017 at 18:37

Subscript[a, 0] = 2;
Subscript[b, 0] = 4;
nN = 1000;
n = 50;
α = 0.05;
n1 = nN*α;
c = 2;


Your use of Table with RandomVariate is unnecessary and inefficient.

at1 = AbsoluteTiming[
SeedRandom[0];
w1orig =
Table[RandomVariate[
BetaDistribution[Subscript[a, 0], Subscript[b, 0]]], {i, nN}, {j, 1,
n}];][[1]];


RandomVariate takes a second argument to generate a List or an Array of variates.

at2 = AbsoluteTiming[
SeedRandom[0];
w1 = RandomVariate[
BetaDistribution[Subscript[a, 0], Subscript[b, 0]], {nN, n}];][[1]];


Comparing the timing

at1/at2

(* 66.6248 *)


Verifying that the random data produced is identical.

w1orig === w1

(* True *)

me = Subscript[a, 0]/(Subscript[a, 0] + Subscript[b, 0]);

se = (Subscript[a, 0] Subscript[b,
0])/((Subscript[a, 0] + Subscript[b, 0])^2 (Subscript[a, 0] +
Subscript[b, 0] + 1));


Note that I have removed the List brackets from the definition of se.

Use of Table to calculate the means and variances is inefficient. Use Map (/@) rather than Table.

at3 = AbsoluteTiming[
m6orig = Table[Mean[w1[[i, All]]], {i, nN}]][[1]];

at4 = AbsoluteTiming[m6 = Mean /@ w1][[1]];

m6orig === m6

(* True *)

at3/at4

(* 24.0071 *)

at5 = AbsoluteTiming[
s6orig = Table[Variance[w1[[i, All]]], {i, nN}]][[1]];

at6 = AbsoluteTiming[s6 = Variance /@ w1][[1]];

s6orig === s6

(* True *)

at5/at6

(* 8.95349 *)

at7 = AbsoluteTiming[e6orig = Table[m6[[i]] - a/(a + b) == 0, {i, nN}]][[1]];


Look at the documentation on Listable and Thread

at8 = AbsoluteTiming[e6 = Thread[m6 - a/(a + b) == 0];][[1]];

e6orig === e6

(* True *)

at7/at8

(* 2.02762 *)

at9 = AbsoluteTiming[e26orig = Table[
s6[[i]] - (a b)/((a + b)^2 (a + b + 1)) == 0,
{i, nN}];][[1]];

s6 - (a b)/((a + b)^2 (a + b + 1)) == 0]][[1]];

e26orig === e26

(* True *)

at9/at10

(* 13.3048 *)

at11 = AbsoluteTiming[
v1orig = Table[
FindRoot[{e6[[i]], e26[[i]]}, {{a, 1.5}, {b, 3.5}},
MaxIterations -> 5], {i, nN}];][[1]];

(* FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 5 iterations.

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 5 iterations.

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 5 iterations.

General::stop: Further output of FindRoot::cvmit will be suppressed during this calculation. *)


However, you can bypass many of the previous calculations (means, variances, equations) and directly calculate the parameters from w1 using FindDistributionParameters

at12 = AbsoluteTiming[
v1 = FindDistributionParameters[#, BetaDistribution[a, b]] & /@ w1;][[1]];


Looking at the first few sets of parameters from each

v1orig[[1 ;; 5]]

(* {{a -> 2.17885, b -> 4.47931}, {a -> 1.60802, b -> 2.75178}, {a -> 1.59492,
b -> 3.24896}, {a -> 2.61362, b -> 5.86455}, {a -> 2.00275, b -> 3.52733}} *)

v1orig[[1 ;; 5]]

(* {{a -> 2.17885, b -> 4.47931}, {a -> 1.60802, b -> 2.75178}, {a -> 1.59492,
b -> 3.24896}, {a -> 2.61362, b -> 5.86455}, {a -> 2.00275, b -> 3.52733}} *)

(at3 + at5 + at7 + at9 + at11)/at12

(* 1.50998 *)


The overall timings are relatively close; however, use of FindDistributionParameters eliminates a lot of unnecessary work.

You should also look at the documentation for FindDistribution and DistributionFitTest