% How many Monte Carlo samples required to achieve a 95% confidence interval
% equal to alpha?  For example, how many Monte Carlo simulations would be
% required to distinguish a p value of 0.001 from 0.002. An application of
% Theorem 5.3.2 on page 331 of
% Larsen & Marx (2001). The application of the decision rule in a Monte Carlo
% trial can be viewed as the outcome of a Bernoulli trial. For most applications
% of the Monte Carlo method, there is an observed pattern. Each Monte Carlo
% simulation is evaluated to determine whether it produces a pattern equal
% to or more extreme than the observed pattern. That probability, which is
% a Bernoulli trial is a binomial process. Theorem 5.3.2 calculates the n
% required to have at least a 100(1-alpha) probability of obtaining an
% observed p value within a distance of d of p.
% Last revised 3/4/11

d=logspace(-6,-1)';% Create 50 logarithmically spaced elements from 10^-6 to 10^-1

% Power=95%, i.e., 95% chance of distinguishing a difference of size d.
alpha=0.05;
n=LMtheorem050302_4th(alpha,d);   % Theorem 5.3.2 from Larsen & Marx, p. 331 in 3rd edition and 373 in 4th
%                              % A result in a Monte Carlo
                               % trial can be regarded as an outcome of a Bernoulli trial
X=[ones(size(d)) log(d)];      % set up the explanatory variable ln(d) for a log-log regression;
                                % A column of all 1's first to fit the Y intercept
Y=log(n);                      % ln(number of trials) from Theorem 5.3.2
B=X\Y                          % Matlab least squares regression

loglog(d,n);xlabel('d in Theorem 5.3.2','FontSize',20);
ylabel('Number of Monte Carlo simulations required','FontSize',20);
title(sprintf('N = %4.2f*{alpha}^{%4.2f}  Power=%4.1f%%',exp(B(1)),B(2),100*(1-alpha)),'FontSize',22);
figure(gcf)
pause

% The previous analysis was for a worst-case scenario, where the p values was 0.5. It would be the
% appropriate equation if you wanted to distinguish a p value of 0.500 from 0.501. But, in testing
% the null hypothesis that 0.005 is different from 0.006, the above
% approach is too conservative

alph=[0.05 0.01 0.001];
for i=1:length(alph)
    alpha=alph(i);
    n=LMtheorem050302_4th(alpha,d,d);   % Theorem 5.3.2 from Larsen & Marx, p. 331. A result in a Monte Carlo
                                   % trial can be regarded as an outcome of a Bernoulli trial
    X=[ones(size(d)) log(d)];      % set up the explanatory variable ln(d) for a log-log regression;
                                    % A column of all 1's first to fit the Y intercept
    Y=log(n);                      % ln(number of trials) from Theorem 5.3.2
    B=X\Y;                          % Matlab least squares regression

    loglog(d,n);xlabel('d in Theorem 5.3.2','FontSize',20);
    ylabel('Number of Monte Carlo simulations required','FontSize',20);
    title(sprintf('N = %4.2f*{alpha}^{%4.2f}  Power=%4.1f%%',exp(B(1)),B(2),100*(1-alpha)),'FontSize',22);
    figure(gcf)
    pause
end
