Run Block of Code in For Loop

Question

Suppose I want to conduct an experiment 10,000 times in which a random sample of size $n=20$ from a beta distribution is collected and then the values of four different estimators are collected. At the end of the 10,000 experiments, I should have a different list with 10,000 instances of each estimator (so essentially I'm trying to get 10,000 random samples for all four estimators). So I started with something like this:

n = 20
mypdf = PDF[BetaDistribution[5, 1], x]
randsamplesize20 = List[]
For[i = 0, i < n, i++, 
 AppendTo[randsamplesize20, 
  RandomVariate[ProbabilityDistribution[mypdf, {x, 0, 1}]]]]
Xbar = Mean[randsamplesize20]
MOM = Xbar/(1-Xbar)
LogSum = Total[Log[randsamplesize20]]
MLE = (-n/LogSum)
Bayes = (n*(n+1))/(1-n*LogSum)
Estimator4 = ((n-1)MLE)/n

But you notice that this code only conducts one experiment. I need to see the average of MOM, MLE, Bayes, and Estimator4 over 10,000 independent experiments. I'm really not very familiar with Mathematica yet, so I tried to wrap all of this in a for loop and ended up with several errors. So my question is, how can I conduct this experiment 10,000 times while still preserving the values of the estimators in each experiment so that I can perform statistical analysis on the random sample of my estimators?

Answer 1

You should use the 2-arg version of RandomVariate. For example, here are 10 samples with size 5:

data = RandomVariate[ProbabilityDistribution[mypdf,{x,0,1}], {10, 5}]

{{0.907865, 0.907415, 0.781752, 0.955632, 0.784422}, {0.941882, 0.931223, 0.929205, 0.621706, 0.945325}, {0.697758, 0.999432, 0.7415, 0.876495, 0.724027}, {0.964263, 0.764248, 0.953871, 0.669235, 0.917145}, {0.654835, 0.665781, 0.839801, 0.880818, 0.992197}, {0.94128, 0.626264, 0.827227, 0.924438, 0.965325}, {0.933157, 0.931406, 0.945465, 0.944415, 0.978303}, {0.870277, 0.785106, 0.95533, 0.974052, 0.95279}, {0.89226, 0.657236, 0.88577, 0.764565, 0.980269}, {0.833877, 0.963077, 0.534098, 0.956507, 0.694219}}

You would use {10000, 20} if you wanted 10000 samples of size 20. Then, you just need to map the appropriate functions over the data. For example, to get your MOM and MLE statistics:

data = RandomVariate[ProbabilityDistribution[mypdf,{x,0,1}], {10000, 20}];

(* MOM *)
Mean[Mean[#]/(1-Mean[#])& /@ data]

(* MLE *)
Mean[-20/Total[Log[#]]& /@ data]

5.21197

5.24855

The other statistics can be determined in a similar way.