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;

psimax = 0;

likelihoodmax = 0;

for B = 1:10000

psi(B) = B/100;

likelihood(B) = 0;

starting = 0;

if (a + c - n) > 0

starting = (a + c - n);

end

ending = m;

if (a + c) < m

ending = (a + c);

end

for A = starting: ending

likelihood(B) = likelihood(B) + nchoosek(m, A)*nchoosek(n, (a + c - A))*psi(B)^(A-a);

end

likelihood(B) = likelihood(B)^(-1);

if likelihood(B) > likelihoodmax

likelihoodmax = likelihood(B);

psimax = psi(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 = psi(B);

end

end

if begin32 == -99

if likelihood(B) > 1/32

begin32 = psi(B);

end

end

if and(begin8 ~= -99, end8 == -99)

if likelihood(B) < 1/8

end8 = psi(B);

end

end

if and(begin32 ~= -99, end32 == -99)

if likelihood(B) < 1/32

end32 = psi(B);

end

end

end

psi1 = input('What is the psi assumed by the first hypothesis? ');

psi2 = input('What is the psi assumed by the second hypothesis? ');

likelihoodratio = likelihood(psi1*100)/likelihood(psi2*100)

psimax

begin32

end32

begin8

end8

plottools