dat = betarnd(3,3,1000,1);
mean(dat)
var(dat)
skewness(dat)

ns = [100,500,1000,5000,10000,50000]
nr = 200;

datTotal = [];
for i=ns
    vr = [];
    rndSamples = betarnd(3,3,i,nr);
    datTotal = [datTotal;mean(rndSamples)]
end

anova1(datTotal')