How can I improve performance of my simulation code?

Question

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[
    RandomVariate[BetaDistribution[Subscript[a, 0], Subscript[b, 0]]], 
    {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[
    RandomVariate[BetaDistribution[OverHat[a][[i]], OverHat[b][[i]]]], 
    {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]]}

Answer 1

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]];

at10 = AbsoluteTiming[e26 = Thread[
      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