%% Hypothesis Testing in MATLAB
%
% Prof. W. R. Ashurst,
% CHEN 3600 FA11,
% Copyright 2011
% Modified: Prof. T. D. Placek
% CHEN 3600 Spring 2012
% Copyright 2012

%% Default Initialization

format long; format compact;
clear all; clc; close all force; tic

%% Tails of the z-distribution
z = linspace(-4, 4);
y = normpdf(z, 0, 1);
plot(z,y)
alpha = 0.05;
a = norminv(alpha/2, 0, 1);
z_patch = linspace(-a, a);
y_patch = normpdf(z_patch, 0, 1);
patch([-a, z_patch, a], [0, y_patch, 0], 'r', 'FaceAlpha', 0.15 )
xlabel('z'); ylabel('P(z)');
a

% The shaded area represents the confidence interval of the z-statistic
% associated with the 95% confidence level.

%% Left tail of the z-distribution
z = linspace(-4, 4);
y = normpdf(z, 0, 1);
plot(z,y)
alpha = 0.05;
a = norminv(alpha, 0, 1);
z_patch = linspace(-4, a);
y_patch = normpdf(z_patch, 0, 1);
patch([-4, z_patch, a], [0, y_patch, 0], 'r', 'FaceAlpha', 0.15 )
xlabel('z'); ylabel('P(z)');
a

% The shaded area represents the confidence interval of the z-statistic
% associated with the 95% confidence level.

%% Right tail of the z-distribution
z = linspace(-4, 4);
y = normpdf(z, 0, 1);
plot(z,y)
alpha = 0.05;
a = norminv(1-alpha, 0, 1);
z_patch = linspace(a, 4);
y_patch = normpdf(z_patch, 0, 1);
patch([a, z_patch, 4], [0, y_patch, 0], 'r', 'FaceAlpha', 0.15 )
xlabel('z'); ylabel('P(z)');
a

% The shaded area represents the confidence interval of the z-statistic
% associated with the 95% confidence level.
%% Sample some data; Repeat 10 times
clc
% Make up 10 data points 'randomly' and plot them on the z-distribution.
alpha = 0.05;
a = norminv(alpha/2, 0, 1);
z_patch = linspace(-a, a);
y_patch = normpdf(z_patch, 0, 1);

for k = 1 : 10
clf
hold on
mu_d = 8*(0.5-rand) % some mean between -4 and 4
sig_d = 0.3+rand; % some standard deviation between 0.3 and 1.3
NN_d = 10;
patch([-a, z_patch, a], [0, y_patch, 0], 'r', 'FaceAlpha', 0.15 )
xlabel('z'); ylabel('P(z)');
axis([-4 4 0 0.4]);
rnd_data=normrnd(mu_d, sig_d, 1, NN_d)
plot(rnd_data, zeros(1,NN_d), 'o',mean(rnd_data), 0, 'x') % x marks mean
pause (2)
end
hold off
%% Tim's calculation of confidence interval
clc
clear x
x = [0.8403, 0.8363, 0.8447];
[MUHAT,SIGMAHAT,MUCI,SIGMACI]=normfit(x, 0.01)
%% Manual calculation of confidence interval
x=[0.8403, 0.8363, 0.8447];
sig=0.0068;
[est1 est2]=CIMeanKnownSigma(x,sig,0.01)
[est3 est4]=CIMeanUnknownSigma(x,0.01)
%% An example z-test for the mean
%{
A laboratory analyzes a sample for an active pharmaceutical ingredient. The
process used to measure the active ingredient concentration has a known
standard deviation of 0.0068 units. The lab measured the same sample of
material three times and found concentration of 0.8403, 0.8363, 0.8447
units. Find the 99% confidence interval for the mean active ingredient
concentration.
%}
clc
clf
ai = [0.8403, 0.8363, 0.8447];
sig_ai = 0.0068;
[~, ~, ci] = ztest(ai, 0, sig_ai, 0.01)
% assume Ho is that mu = 0, just to run the test and get the confidence
% interval information.

h = ztest(ai, ci(1), sig_ai, 0.01)
h = ztest(ai, ci(2), sig_ai, 0.01)

%{
The confidence interval is found to be [0.83032, 0.85055] units. Thus, any
hypothesis test where the null hypothesis involves a mean inside this
confidence interval will fail to reject (FTR) the null.
%}

[h, ~, ~, z_star] = ztest(ai, mean(ci), sig_ai, 0.01)

%% Another Example z-test for the mean

%{
Consider a simple random sample of 200 exectutives whose mean blood
pressure is 126.07 mmHg. Assume that the "normal" population mean is 128
mmHg and the population standard deviation is 15 mmHg. Does the blood
pressure data for the executives differ from the population at the 95%
confidence level? Is it less than the general population at the 95%
confidence level? More than?
%}

clc; clear all;
NN = 200;
bp = 126.07*ones(1,NN);
mu_bp = 128; sig_bp = 15;
alpha = 0.05;

[h, p, ci, z_star] = ztest(bp, mu_bp, sig_bp, alpha, 'Both')
% FTR Null (h == 0) - The data is not strong enough to reject the idea
% that there is no difference in the blood pressure for the exectuives
% versus the normal population.

[h, p, ci, z_star] = ztest(bp, mu_bp, sig_bp, alpha, 'Left')
% Reject the Null (h == 1) - The data is strong enough to reject the
% assertion that the blood pressure for executives is no different than the
% normal population in view of the alternative that is it less than the
% normal population. You beleive, then, that the blood pressure for
% executives is "on average" lower than the blood pressure for the normal
% population.

[h, p, ci, z_star] = ztest(bp, mu_bp, sig_bp, alpha, 'Right')
% FTR Null (h == 0) - The data is certainly not strong enough to reject the
% idea that the blood pressure for executives is higher than the normal
% population.

%% Executive blood pressure - the "Hard Way"
%
% This time, we manually compute the z-statistic and p-value and compare
% the p-value to alpha, according to the "tailedness" of the test we are
% conducting.

z_star = (mean(bp)-mu_bp)/(sig_bp/sqrt(length(bp)))

% Ha: mu1 ~= mu2
two_p = 2*normcdf(z_star, 0, 1)
% Since 2P > alpha, we Fail To Reject the Null.

% Ha: mu1 < mu2
p_left = normcdf(z_star, 0, 1)
% since p_left < alpha, we reject the null.

% Ha: mu1 > mu2
p_right = 1 - normcdf(z_star,0,1)
% since p_right > alpha, we Fail To Reject the null.

%% Less Executives

%{
Now consider a simple random sample of only 100 exectutives whose mean
blood pressure is 126.07 mmHg. As before, assume that the "normal"
population mean is 128 mmHg and the population standard deviation is 15
mmHg. Does the blood pressure data for the executives differ from the
population at the 95% confidence level? Is it less than the general
population at the 95% confidence level? More than?
%}

clc; clear all;
NN = 100;
bp = 126.07*ones(1,NN);
mu_bp = 128; sig_bp = 15;
alpha = 0.05;

[h, p, ci, z_star] = ztest(bp, mu_bp, sig_bp, alpha, 'Both')
% FTR Null

[h, p, ci, z_star] = ztest(bp, mu_bp, sig_bp, alpha, 'Left')
% FTR Null

[h, p, ci, z_star] = ztest(bp, mu_bp, sig_bp, alpha, 'Right')
% FTR Null

% This time, with only 100 executives in the measured pool, the data is not
% strong enough to reject the null in view of any of the alternative
% hypotheses!

%% The Speed of Light - an example t-test
%{
Around 1880, Simon Newcomb attempted to estimate the speed of light by
measuring how long it took light to pass from his laboratory to a mirror
and back, a distance of almost 7400 meters. Newcomb's data is represented
in the vector below. Modern measurements of the speed of light would
translate to a measurement of 33.02 in Newcomb's experiment. Is the modern
estimate significantly different at the 99% confidence level?
%}

clear all; clc;
Newcomb = [28, 22, 36, 26, 28, 28, 26, 24, 32, 30, 27, 24, ...
33, 21, 36, 32, 31, 25, 24, 25, 28, 36, 27, 32, ...
34, 30, 25, 26, 26, 25, 23, 21, 30, 33, 29, ...
27, 29, 28, 22, 26, 27, 16, 31, 29, 36, 32, 28, ...
40, 19, 37, 23, 32, 29, 24, 25, 27, 24, 16, ...
29, 20, 28, 27, 39, 23];
% Newcomb = Newcomb + 24800;
% Newcomb = Newcomb*1e-9;
alpha = 0.01;

[h, p, ci, t_data] = ttest(Newcomb, 0, alpha)
% done to establish the confidence interval for the Newcomb data.

[h, p, ci, t_data] = ttest(Newcomb, 33.02, alpha)
% We reject the null (h == 1) and can say that Newcomb's data is
% significantly different from the modern estimate. We may be wrong about
% this only 1% of the time as it is the Type I error. This is also obvious
% in view of the confidence interval established immediately beforehand.

%% Can Third Graders Read? - An example two sample t-test.

%{
The degree of reading power (DRP) scores for a control group and a treated
group of third grade elementary school pupils is given below. Does the
treated group show an improvement in their DRP at the 95% confidence level?
What about the 99% confidence level?
%}

clear all; clc

cd = [42, 46, 43, 10, 55, 17, 26, 60, 62, 53, 37, 42, 33, 37, 41, 42, ...
19, 55, 54, 28, 20, 48, 85];
td = [24, 56, 43, 59, 58, 52, 71, 62, 43, 54, 49, 57, 61, 33, 44, 46, ...
67, 43, 49, 57, 53];

alpha = 0.05;
[h, p, ci, tdata] = ttest2(cd, td, alpha, 'Left', 'Unequal')
% We Reject the Null (h == 1) and claim that the data is strong enough to
% refute the hypothesis of no difference and can assert that the treated
% group shows a difference in DRP at the 95% confidence level.

alpha = 0.01;
[h, p, ci, tdata] = ttest2(cd, td, alpha, 'Left', 'Unequal')
% Uh-oh - FTR NUll (h == 0). The data is not strong enough to refute the
% hypothesis of no difference at the 99% confidence level.

%% Two Sample t-test without the measured data!

%{
Many times you may be provided with "digest" information about a data set.
For example there were N samples, the mean was X and the standard deviation
was s, and the data was normally distributed. This can be convenient when
there are hundreds of data points. We may not use the MATLAB functions
such as 'ttest' and 'ttest2' to deal with these scenarios. We can still
conduct the hypothesis test the "hard way" by manually computing the
challenge statistic and p-value and manually comparing the p-value to
alpha.

Assume that you have a sample of response digests from Males and Females
about some social insight test. The relevant data are loaded in the script
immediately below this comment.
%}

clear all; clc
M_n = 133; M_mu = 25.34; M_sig = 5.05;
F_n = 162; F_mu = 24.94; F_sig = 5.44;

% Do the Males differ in opinion from the Females at the 95% confidence
% level?

alpha = 0.05;
t_star = (F_mu - M_mu)/(sqrt(F_sig^2/F_n + M_sig^2/M_n))
p = 2*tcdf(t_star, min(F_n, M_n)-1)

% Since p > alpha, we fail to reject the null. Note that since this is a
% two-tailed test, we double the native tcdf value. Since the p-value is
% so high, the data is not nearly strong enough to refute the hypothesis of
% no difference between the samples.

%% Two Sample t-test without the measured data - Using Monte Carlo (a.k.a. the "really hard way")

%{
Repeating the example from above, but instead of using the formulas for the
t-test we will construct an artifical data set using Monte Carlo methods.
We will have to generate the data set and then interrogate that data set to
see if its characteristics are close enough to the provided information.
For the acceptable data sets we then carry out the normal hypothesis tests,
keeping track of how many times we reject the null. As usual in Monte
Carlo methods, we need to conduct the simulation many times to establish
the "long-term" behavior. Note that mean_tol and sig_tol are arbitrary
guesses and can be refined as needed. The lower these values are, however,
the more times data sets will have to be generated before "acceptable" data
sets are discovered.
%}

count_valid_simulations = 0;
count_simulations = 0;
count_reject_null = 0;
mean_tol = 0.5; sig_tol = 0.05;

while (count_valid_simulations < 100 && count_simulations < 5e6)
M_data = normrnd(M_mu, M_sig, 1, M_n);
F_data = normrnd(F_mu, F_sig, 1, F_n);
if (abs(mean(M_data) - M_mu) < mean_tol && ...
abs(mean(F_data) - F_mu) < mean_tol && ...
abs(std(M_data) - M_sig) < sig_tol && ...
abs(std(F_data) - F_sig) < sig_tol)
count_valid_simulations = count_valid_simulations + 1;
[h, p, ci, tdata] = ttest2(M_data, F_data, alpha, ...
'Both', 'Unequal');
if (int8(h))
count_reject_null = count_reject_null + 1;
end
end
count_simulations = count_simulations + 1;
end
fraction_of_reject_nulls = count_reject_null/count_valid_simulations
M_pl = hist(M_data, 20); F_pl = hist(F_data, 20);
bar([M_pl', F_pl']);
legend('Males','Females')
%{
Since the 'count_reject_null' is typically 2 or less out of 100 (with these
modest _tol values), we may say that the simulated data pools that
correspond with the characteristics provided consistently fail to allow us
to reject the null. In other words, the data characteristics are not
strong enough evidence that the Males differ in opinion from the Females.
This is especially true if we set the 95% confidence interval, since
'count_reject_null' has not been observed to be greater than 3 out of 100
tries, and we are allowed to be wrong 5 out of 100 tries. Of course, if we
restrict the mu_tol and sig_tol further, the observations of such "reject
nulls" *should* be less, in view of our analytic test above. This fact is
borne out by experimentation, but at the expense of much longer computation
times.
%}
toc