SEd = input('What is the sample standard error? ');
meand = input('What is the sample mean? ');
n = input('How many subjects were run? ');

Vard = n*SEd^2;
SSd = Vard*(n-1);

likelihoodmax = 0;
theta(1) = meand - 5*SEd;
inc = SEd/100;

for B = 1:1000
theta(B) = theta(1) + (B-1)*inc;
likelihood(B) = (SSd + n*(meand - theta(B))^2)^(-n/2);

if likelihood(B) > likelihoodmax
likelihoodmax = likelihood(B);
end
end
for B = 1:1000
likelihood(B) = likelihood(B)/likelihoodmax;
end

outofrange = meand - 6*SEd;

begin8 = outofrange;
begin32 = outofrange;
end8 = outofrange;
end32 = outofrange;
for B = 1:1000
if begin8 == outofrange
if likelihood(B) > 1/8
begin8 = theta(B);
end
end
if begin32 == outofrange
if likelihood(B) > 1/32
begin32 = theta(B);
end
end
if and(begin8 ~= outofrange, end8 == outofrange)
if likelihood(B) < 1/8
end8 = theta(B);
end
end
if and(begin32 ~= outofrange, end32 == outofrange)
if likelihood(B) < 1/32
end32 = theta(B);
end
end
end

theta1 = input('What is the population mean assumed by the first hypothesis? ');
theta2 = input('What is the population mean assumed by the second hypothesis? ');
B1 = int16((theta1 - theta(1))/inc  + 1);
B2 = int16((theta2 - theta(1))/inc  + 1);
likelihoodratio  = likelihood(B1)/likelihood(B2)

begin32
end32
begin8
end8