a = input('What is a? ');
b = input('What is b? ');
c = input('What is c? ');
d = input('What is d? ');
 
m = a + b;
n = c + d;
k = a + c;
l = b + d;
 
gammamax = 0;
likelihoodmax = 0;
 
 for B = 1:10000
    gamma(B) = B/100;
    likelihood(B) = 0;
    for A = b:(m+n-d)
        likelihood(B) = likelihood(B) + nchoosek(A-1, b-1)*nchoosek(m+n-A-1, d-1)*gamma(B)^(A-m);
    end
    likelihood(B) = likelihood(B)^(-1);
    if likelihood(B) > likelihoodmax
        likelihoodmax = likelihood(B);
        gammamax = gamma(B);
    end
 end
 for B = 1:10000
     likelihood(B) = likelihood(B)/likelihoodmax;
 end
 begin8 = -99;
 begin32 = -99;
 end8 = -99;
 end32 = -99;
 for B = 1:10000
     if begin8 == -99
         if likelihood(B) > 1/8
             begin8 = gamma(B);
         end
     end
     if begin32 == -99
         if likelihood(B) > 1/32
             begin32 = gamma(B);
         end
     end
     if and(begin8 ~= -99, end8 == -99)
         if likelihood(B) < 1/8
             end8 = gamma(B);
         end
     end
     if and(begin32 ~= -99, end32 == -99)
         if likelihood(B) < 1/32
             end32 = gamma(B);
         end
     end
 end
 
 
  gamma1 = input('What is the gamma assumed by the first hypothesis? ');
  gamma2 = input('What is the  gamma assumed by the second hypothesis? ');
  likelihoodratio  = likelihood(gamma1*100)/likelihood(gamma2*100)
    
     gammamax
     begin32
     end32
     begin8
     end8