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