0
$\begingroup$

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]]}
$\endgroup$
  • 5
    $\begingroup$ 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/… $\endgroup$ – gogoolplex Oct 15 '17 at 13:16
  • $\begingroup$ in my simulation step 5 takes the most time,precisely when calculate SuperStar[T] $\endgroup$ – ola farouk Oct 15 '17 at 18:37
6
$\begingroup$
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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.