function nag_stats_sol()

% Compute chi-sqaured statistic wuth G01AF
nobs_d = [0 2 2 3 2; 0 6 1 7 0; 3 2 1 3 0];
nobs = nag_int(nobs_d)

[nobsOut, num, pred, chis, p, npos, ndf, m1, n1, ifail] = ...
                     g01af(nobs);

tail = 'U';
df = double(ndf);
[p, ifail] = g01ec(chis, df, 'tail', tail);

fprintf('Chi-squared statistic using G01AF gives p = %g \n',p)
test_p(p)

nsim = 2000;
rows = nag_int( sum(nobs,2) );
cols = nag_int( sum(nobs,1) );

p = monte_chi(chis, rows, cols, nsim);
fprintf('Chi-squared statistic using Monte Carlo gives p = %g \n',p)
test_p(p)


function p = monte_chi(ochis, rows, cols, nsim)
% Compute chi-sqaured statistic with Monte Carlo method 
   mode = nag_int(2);
   igen = nag_int(1);
   iseed = [nag_int(1762543);
            nag_int(9324783);
            nag_int(42344);
            nag_int(742355)];
   r = zeros(100, 1);
   
   cnt = 0;
   
   for i = 1:nsim
     [nobs, iseed, r, ifail] = g05qd(mode, rows, cols, igen, iseed, r);
     [nobs, num, pred, chis, p, npos, ndf, m1, n1, ifail] = ...
                     g01af(nobs);

     if chis > ochis
         cnt = cnt + 1;
     end
   end
   
   p = (cnt+1) / (nsim+1);

end

function test_p(p)
% Test the null hypothesis
   if p < 0.05 
       fprintf('Reject null hypothesis, there is an association\n\n')
   else
       fprintf('Accept null hypothesis, there is no association\n\n')
   end
end

end