u=rand(1,1000000);
est=u.^(1/3).*(1-u).^(1/4);
m=mean(est);
se=sqrt(var(est)/length(est));
disp(['mean of 100000 samples ' num2str(m) '  standard error ' num2str(se)])
%***************************************
disp(' press any key to continue: ESTIMATING PI')
pause
u=rand(1,100000);
u=rand(1000000,2);
est=u(:,1).^2+u(:,2).^2 <1;
disp(['estimate of pi is ' num2str(4*mean(est))])
se=4*sqrt(var(est)/length(est))