function observed=binomialBayes(N,p,alpha,beta,observed)


%H0: p=p0... Note: Null distribution is using the true p (in practice, you only wish you could do this)!
%Ha: p<p0


for(i=1:N+1)
    bpdf(i)=binopdf(i-1,N,p);
end

%observed: if NaN, generating a value.... 
%then you can reset alpha and beta to show the prior impact
if(isnan(observed))
    observed=binornd(N,p);
end

clf
subplot(1,3,1)
bar(0:N,bpdf)
ylabel('Probability Mass')
xlabel('# Of Observed Heads')
hold on

bar(observed,binopdf(observed,N,p),'r')

pVal=binocdf(observed,N,p)  %changed from old alternaltive which corresponded to H_a p>p_0... now set for p<=p_0 with :binocdf(observed,N,p)
title(['p-value = ',num2str(pVal)])

vals=100;
for(i=1:vals+1)
    paramVal(i) = 1/vals*(i-1);
    like(i)=paramVal(i)^observed*(1-paramVal(i))^(N-observed);
end
like=like*vals/sum(like);
subplot(1,3,2)
hold on
plot(paramVal,like,'r','linewidth',2)
plot([p,p],[0,3],'b','linewidth',2)
ylabel('Likelihood')
xlabel('p')
legend('Likelihood (not a density!)','Truth')

for(i=1:vals+1)
    paramVal(i) = 1/vals*(i-1);
    post(i)=betapdf(paramVal(i),alpha+observed,beta+N-observed);
    prior(i)=betapdf(paramVal(i),alpha,beta);
end

subplot(1,3,3)
axis([0,1,0,10])
hold on
xlabel('p')
plot(paramVal(2:vals),post(2:vals)*vals/sum(post(2:vals)),'r','linewidth',2)
plot(paramVal(2:vals),prior(2:vals)*vals/sum(prior(2:vals)),'k--','linewidth',1)
legend('Posterior Density','Prior Density')
