
% Binomial Estimation and Testing: Proportions from the Same Sample

%A 8196-long nucleotide sequence from
%spider monkey Ateles geoffroyi,
 %obtained from EMLB Nucleotide sequence alignment DNA database
 %
dna = ['gatctcttcttcgtggttagtattctttgtgttctgctaaagaactctttgcctaccctc'... 
 'aggtttgcttaggtctacaatctactgtgagttgatttttttaatggacagtgagatgta'...
 'ggaattgtttctttcttatggatttccagtcattcagcaccatttattaaaaagagcccc'...
 'ctgtgagcccctggtcttggtgccttctgtgttataaatcaggtgttcctgtgtgtctat'...
 'gcgtgggcccatatgtatatttctggactctcttttgccccactggtctctttatctatc'...
 'caagtgataacccctaccaccatcactccttcagctttataatgtcttgatagtcaatag'...
 'agtatatatttcaactttgttcttcaagagtgtcttgataattctcagtcctttgcattt'...
 'tcacataaaacttaaaatctgccttccatttgcatacacgcacacacacaaaaagaaaac'...
 'tctgctagaaatttgatagggatcactttgaatctctctatcaatttggggagactttgc'...
 'aatgctgagtcttccaatccatgaacatgatatatgcttccatttatgcttcaatttttc'...
 'ttaataacttttaaagctttcagtgtcaagttcttatacatattttattagataaattac'...
 'aggtgtttggttttctactataaatagaatcattttgaaattatcttatttgttgctgat'...
 'atatggaaatactgactttgtaaatatgcagtcttactcaatgcacttattagttctggt'...
 'tatttatttgtagaataaacaaaatcatatcatttgtaaaaaaacaggggggttgggttt'...
 'tgtttgtttgtttgtttgtttgagacagggtctcactctgtcatgcagactggagtggag'...
 'tggcacagtcatggctcactgcagccttgacctcttgggctcaagcctcctgagtagctg'...
 'cgaccacaggaatacaccatgaccaagccttgcctgatttttttttttttcccagaaaca'...
 'gggtttctctatgttggccaggctggtcttgaactcctgtgctcaagcagtcctcccgcc'...
 'ctggcctcccaaagtgctgggattacaggcatgagccaccgcacctggtccagttttatt'...
 'tattctctttcagttcttactggctttatttcatgtttgcttttactgaactggacagga'...
 'ttttgagcataatgctgatgtcataatagtagacattaaaagaaaagcacgtatattctt'...
 'accactaaattttatgtttggccaggcatggtggctcatgcctgtaatcccagcactttg'...
 'ggaggccgaggcaggtggatcacctgaggtcaggagcttgagaccagcctggccaacatg'...
 'gtgaaaccctgtctctactaaaaatacaaaaactagctaggtgtgttagcacatgcctgt'...
 'aatctcagctactagggaggctgaggcagcagaattgcttgaacccaggggtcagaggtt'...
 'gcagtgacccaagattgcactgctgcacttgaacctgggcaacaaagcgagactctgtct'...
 'caaaaataaatttaaaaaaaagttatgtttgctatggagttgtttttgttgcttctgctg'...
 'cttcatacatttcctttatcaggttatagaattttccttctatttttctttttatctttt'...
 'ttttttttttgagacaaagtcttgctctgttgcccaggctggagcgcagtggagtgatct'...
 'cagctcaccataacctctgcctcccaggttcaagcgattctcatgcctcagcctcccaag'...
 'tagctgggactacagcatgcgccaccacacccagctaatttttgcatttttagtagacac'...
 'agggtttcaccatgttggccagcatggccttgaactcctggcctcaggtgatctgcccgc'...
 'ctcgacctcccaaagtgttgggattacaggcgtgagccactgtggctgaccccttctatt'...
 'tttatttgccaagagtttttattatgaatagatggcagattttatcaaatgctttttctg'...
 'catctgttgaaataatcatgtgactttcctcctttattcagttattgtagtgaattaaaa'...
 'ttattttattttcaaatgttacaccaaacttgtattccttgaataaatctaattcagctg'...
 'tgatatattagcattttatacaagtatggctaggttctatttgctaatattttgcttaga'...
 'acttttgtatctatgacataaacgagactggcttattatttactctgcttataatatcct'...
 'tattttggtatcaaggttattcttggcctcataaatcaagttgaaaagcattttctcatc'...
 'tctggaaagttttaataagattggcagttttccttaaatgtcagatagaatccactggtg'...
 'aagccaaatgggtctggaaatttctttattggaaagtttcaaattaaatattccatttat'...
 'ttaacctagaagactattcatattttctatttctttttttgtcatttttgctaagttgtg'...
 'tttttctaggaacttgcctgttatacttaggctttcgcatatatctgcatgatcgcttgg'...
 'gtttttacagctgtaggctctgtaataatgtcacccaccctttcctgattctgattacct'...
 'gcacttctctcttattccttgaacattcttcaacaggggctgggtttcttagattcccag'...
 'gaactgtttctccctctctctctctccctctttctccccacttcctctatctgcagggat'...
 'ttaactggactctgttcatacaaactgtgctatcctctgtcaaaatcaagctgctgccag'...
 'atgaggaagtggtggtctatggcatcccctacctgcagaaccttgaaaacatcatcgaca'...
 'cctactcagccaggtgaggcccagggagggtctccaaaggcctcagcatctttgcccctg'...
 'ctctcccatcctccccaacacagcttctgctcccaggaggggaccctatgaagcatctgc'...
 'cctgaggggtgtgggacatgaggggttaagctttagtggccagtggagtcttacattaaa'...
 'aatgtgtcacatcccaagcacccggggcagatgggctggatcttggaggtggttgtcctg'...
 'gaggttttctccagcacccataactctctggggtggaacgtgctggccaagggccctaag'...
 'ccatatctcctcccttgcaaagacctctcacctgttctcccctcaccctctaccttctgc'...
 'ccccaggaccatacagaactacctggtctggcgcctggtgctggaccgcattggtagcct'...
 'aagccagagattcaaggacacacgagtgaactaccgcaaggtgagccccttccctgccca'...
 'gagcctggccgggcactcattcattcaggccacgatgctgagcataagctgtaggcaaac'...
 'accgaggtccacaccagggtcacagccatggacgaggcagggcgaggcgcccaccccagt'...
 'ggggagtggaccatgacagtcaacaaggcagaacgtagctctgccaaaaccgcagagtat'...
 'gatgtggcggggaggcagggagggacgaggacagggagcactgctccggggaggcgacat'...
 'ttgagaccttaaaggcctgaagaggcggcggggtgatttcccgggagcagaaacagccag'...
 'tgccagtgcagaggccagaggcaggagcaggcaggggccggacccgtggtcagcagcgga'...
 'caggaaggagaagatggcagggagggaggcgaggtgaccacgggccacagcagccggttt'...
 'ggattgtactccctgtctgaaagccgcaggggtgctaagagcagagacgcggtcagatat'...
 'gctgctcctgtagtcactgtggctgctccagggctggggtgagggtggagacaaggaacc'...
 'gcagggaggccctagggggattctggctgttctggggagagtgagttgagtagggggcgg'...
 'ggtcagcatggaggccgatgtggtcacccaggtgaggtgacaggcttggccatggtgctg'...
 'gcagacacacagacgcaggggaggtgttgaccgcagcagctgagggattgggtatagggg'...
 'ctaagggcagaggggagtcgggactctcaagacctttgttagagtaatgggagggacaga'...
 'ggtccccactgcaatggggagactgaggcagaggcaagctgggggcaagtcaagcctcct'...
 'gttgccatgggggggtaggtcagagacatctgtgagtcaagcgtcaggggggcactgagt'...
 'ctaggagctttggggaggaagaaagggctgcctccctttcctaaagggcaggaagtcctg'...
 'cagcccagaaggtgtggcaggacgctgggcagagtgctccctgggctttgcctgagtctg'...
 'gggcccagagtctcgtgctcctgtgtccctggggacagatgagcaggtgttctcacagct'...
 'gcttcaacaatggagagaccagggccctgggagctggtggctgagagggacgtggtccca'...
 'ggccttgctctggggtcagcgccagccctgtgccttgagactggggctggagtcctggag'...
 'agcagagtccagccagacaaccaggacaggggccggggaggaacccccccaactcctcag'...
 'ggtaaaggcccagctctcagtctgcggaagttccccaatcacaataccatgcccatcctt'...
 'ccccccgccgcagcagcacccagaggtgaggctgagatcatcccgtgcccccgctagccc'...
 'cggggcaccctgaggctgaggataggtgggtctgctcaccccgtaagcccagggccccac'...
 'gccaggcctgggaaggctgacgtgggggcccgtgggatgatcgggtggccccaccaggcg'...
 'ctgtttggcacaatggtggaggaggtgcgctggcgtgaatgtgtgggctacgtcaacagc'...
 'aacatggagaacgccgtgggctccctctacgtcagggaggcgttccctggagacagcaag'...
 'agcatggtgggcgcccctcatccatggcccacgcggggccagccttccacagggacaccg'...
 'agcctccctgccctgagactcccaccccgcccaggctcaggggcacctgccccagctcac'...
 'cagacagagggaggctaagggaagccgtgcccaccccatgcccctgaccccagccggcca'...
 'ctgggccccagcctgagccctgcaggtccagctgggaaaggcacaggtggtacaagttgc'...
 'ccctggcccacgaggacccacagcccagcagggaggaccctgaagttgggttgtcatctc'...
 'atccccccaccggccccacacctcacctgagaggcccacagctgggcagccctgggcctt'...
 'gccccggcggggctcagaccagacctccccaaagcccttgaggctggtcattcgccctcc'...
 'agatgagtctgggacacgagggccatgcccgattcccaggtcagagaactcattgacaag'...
 'gtgcggacagtgtttgtggagacgctggacgagctgggctggatggacgaggagtccaag'...
 'aagaaggcgcaggagaaggtatgtggctcgcagaggtcagggcagccccgactggatctg'...
 'ggcacgggcagcaccgccccgacaaagcctcctgacggtccccacggggagaaggcagcc'...
 'ccgggaggcgcagacaggaccccacggggagaaggcagccccgggaggcgcagacaggac'...
 'cccacggggagaaggcagccccgggaggcgcagacaggaccactcgcccactgtggagct'...
 'gaggaaaagcaggcaggggcaggaaccaggcctgggacatccccagcaccccccatgagc'...
 'cagcgggggtagcggtggcctcgggagcaagggtcagacccaccagaagctcctctggac'...
 'ccaagtgacagacggcagccgagtgacccctccagggcccctgcccccacccggacccct'...
 'ggcgggctgttgggctctaggtcaggtgggcggagccttccatattctgtccccattaaa'...
 'cctccccagagcctgcaaggtgggggccagatcccggagccccccagcccctcaccacag'...
 'gacagggagagcacagccaggatgccccgggctgacccctcccacagtgtggcactcact'...
 'gaccggcccccgggctgactccatccgcagtggggcacccactgaccagtccctcctccc'...
 'cactcaggccatgagcatccgggagcagatcgggcaccctgactacatcctggaggagat'...
 'gaacaggcgcctggacgaggagtactccaatgtgcgcccccaaccccagccctaaactgc'...
 'cccccacaccgtggacatgggccctgattgcctgggggctcctctgtgccaagctctgtg'...
 'ccctcctcccaggcgctaacacagctgtcctcaccgtgccctagggtggtggctccagca'...
 'ctaccctcagtgtacagatggggaaactgaggcctggagggacgagggacttgctcagag'...
 'tcatgtaaatgcccacagccagggctgccccgttcacctgcacccacagccgactgtcca'...
 'cgtgcggacgcctggcaaggcttgaggtccttccccaaccacccaggacctctccccttg'...
 'ctcctctgagggcccggccccagcatctggcagtgactgcccctctcccgtagctgaact'...
 'tctcagaggacctgtactttgagaacagtctgcagaacctcaaggtgggcgcccagcgga'...
 'gcctcaggaagcttcgggaaaaggtggacccaaatctgtgagtggcatcgctggtgcggc'...
 'atgggcggggggggccagggagtggccagtaaggcgtgcagggctgcctgccaggcgggg'...
 'gactggatgcgtggctggggtccagtgccagctgaggggtgcagtcactgggcctatctt'...
 'ggacagagtcccccatcaccctctgctttctggggcccaagtcgcctcctcctcagggac'...
 'cccaggaaggaaccagatcttgggataggcctgggcccttgctctggtgtggccagggta'...
 'ggccctggtccatcaaggcagagagagcccagagggagggtccccggagaggcggctccc'...
 'agagaccggcttgtctatcacagctggatcatcggggcggcggtggtcaatgcgttctac'...
 'tccccaaaccgaaaccagattggtaagttcctccccctccctgcctgcccaggccctttg'...
 'gtgccagagttgaaggggtggcttgcagccccctgcccatcccttcctcctgggggcacc'...
 'ccagcctgttcttggaacattctcctatgtcaccccttgctgcccccgtggggaacactc'...
 'actcccgtgcccatgagcggcttcgggagggacttggcaggcggcccagagcctttcctg'...
 'caggcagaacaggatggtgggaggaggccggactgaggtcctgggtctgccaggcacgag'...
 'cccgcccagcctctgaacctccgtgtcttcatctgtaaattgagggaatcctgacctctg'...
 'cctcccggtctgtgggatcaagtgtcagaacccctggcagcgccatccctgagcctggca'...
 'cacaggagatgtgtccactcacagcttattctttgtcccgcccttggggacagtattccc'...
 'tgccgggatcctccagccccccttcttcagcaaggagcagccacaggccttgaactttgg'...
 'aggcattgggatggtgatcgggcacgagatcacgcacggctttgacgacaatggtgagcg'...
 'gccgcccctccaatcacacacgccctcccctctccgcggagaagcagccaggaaggggcc'...
 'tgggggcgctccattcaaccccttcacggtttcagtggggacaccgaggccaggtcctca'...
 'gggagctgtgggggcagagcaggctttgaatgccaggctggggccctgtctgggggcact'...
 'gactcccagcatggaggaggggacccagggactccagcccatgaccctgactccacctgg'...
 'cactgggcacttccctcccacggcctgtgacatggggacgacaacactgcccaatacagc'...
 'agcgaagcccagatcccagcctgggagactctccagacggccagcatcctccaccggggg'...
 'gtgtgctacccccgccatggccatgctgacctccgcccctcggtggggactgccctgcct'...
 'gggacaggggccccagctggttttcccctgtg']; 

% (a) Demonstrate that the proportion of 'a' is significantly
% smaller than the proportion of 't'.
% (b) Demonstrate that the proportions of 'g' and 'c' are not
% significantly different.

%
%-----------------------------------------------
% Point estimates
n = length(dna);                       %8192
pahat = sum( dna == 'a')/n; %0.2144
pghat = sum( dna == 'g')/n; %0.2769
pchat = sum( dna == 'c')/n;  %0.2787
pthat = sum( dna == 't')/n;   %0.2301


% (a) H0: pa = pt   H1: pa < pt
% (b) H0: pg = pc   H1: pg ~= pc

% Var(pa - pt) = pa*qa/n + pt*qt/n + 2*pa*pt/n
%              = [pa + pt - (pa - pt)^2]/n
%              = (pa + pt)/n  [under (a) H0]

%(a)
Zat = sqrt(n) * (pahat - pthat)/sqrt(pahat + pthat); %-2.1379
pat = normcdf(Zat); %0.0163
%(b)
Zgc = sqrt(n) * (pghat - pchat)/sqrt(pghat + pchat); %-0.2224
pgc = 2*normcdf(-abs(Zgc)); %0.8240
% %-----------------------------------------------
  sprintf(' pat= %7.4f', pat)
  sprintf(' pgc= %7.4f', pgc)
% 
% % --------------------------------
% For confidence intervals we use
% std(pahat - pthat)=sqrt{ [pa + pt - (pa - pt)^2]/n },
% evaluated at pa=pahat and pt=pthat.
sat=sqrt((pahat + pthat - (pahat - pthat)^2)/n) %0.0074
[pahat - pthat - sat*norminv(0.975)   ...
       pahat - pthat + sat*norminv(0.975)]
   % -0.0302   -0.0013
sgc=sqrt((pghat + pchat - (pghat - pchat)^2)/n) %0.0082
[pghat - pchat - sgc*norminv(0.975)   ...
       pghat - pchat + sgc*norminv(0.975)]   
    %-0.0180    0.0143
%%
sati=sqrt( (pahat * (1-pahat) + pthat * (1-pthat))/n) % 
[pahat - pthat - sati*norminv(0.975)   ...
       pahat - pthat + sati*norminv(0.975)]
  %
%%
   atc=1756; tac=1885; at=0; actc=2268+2283;
   N=8192;
   A=at ; B= atc; C= tac; D=actc;  %table
   
   p1hat =(A+B)/N; p2hat = (A+C)/N; p11hat=A/N; p12hat=B/N;
   p21hat=C/N; p22hat=D/N;
   
   %the two definitions of sample variance for risk difference
   % coincide
   %(1) by paired tables 
  s2= p1hat * (1-p1hat)/N + p2hat * (1-p2hat)/N  - ...
      2* (p11hat*p22hat - p12hat * p21hat)/N
 
  %(2) by dependent proportions
 u2= (C/N + B/N - (C/N - B/N)^2)/N

 