%tpower.m
%Perform t-test and plot results of Observed Power analysis

%create two data samples
sample1=[6 5 7 9 9 8 6 7 5 9];
sample2=[5 6 5 7 6 4 6 7 5 5];

%conduct one-tailed (sample1 > sample2) t-test (alpha=0.05)
alpha = 0.05;
[h,p,ci,stats]=ttest2(sample1,sample2,alpha,'right');
display(['Null rejected? ',num2str(h)]);
display(['p-level:   ',num2str(p)]);
display(['df:        ',num2str(stats.df)]);
display(['t-stat:    ',num2str(stats.tstat)]);
display(['pooled sd: ',num2str(stats.sd)]);

%groups size and related parameters
n1 = length(sample1);
n2 = length(sample2);
df = n1+n2-2;                               %degrees of freedom
d = (mean(sample1)-mean(sample2))/stats.sd; %Cohen's d (Effect Size)
ncp = d * sqrt((n1*n2)/(n1+n2));            %noncentrality parameter

%critical t-value 
tcrit = tinv(1-alpha, df);
%derive beta
beta = nctcdf(tcrit, df, ncp);
%compute power
ePower = 1-beta;

%plot central and noncentral t distributions
x = [-5:0.1:8];
nct = nctpdf(x,df,ncp);
tdist = tpdf(x,df);
figure, clf, plot(x,nct,'b-','LineWidth',2);
xlabel('t-value'), ylabel('Probability Density');
title('t-Test Power Analysis');
hold on;
plot(x,tdist,'g-','LineWidth',2);
%annotate figure
legend('nct','t');
line([tcrit tcrit],[0 0.5],'LineWidth',2,'Color',[0 0 0]);
fstr(1)={['t-critical=',num2str(tcrit),' (alpha=',num2str(alpha),')']};
fstr(2)={['df=',num2str(df),',  nc-parameter=',num2str(ncp)]};
fstr(3)={['Effect Size (d) = ',num2str(d)]};
fstr(4)={['Beta = ',num2str(beta)]};
fstr(5)={['Observed Power = ',num2str(ePower)]};
text(-5.75,0.4,fstr);
hold off;
