function [t0, t1, p, n, x] = BIET(delta, c, alpha, beta, traces_res)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIET - Function of bayesian interval estimation testing
%
%   BIET(DELTA, C, ALPHA, BETA, TRACES_RES)
%
%   DELTA - half size of an estimation interval which may contain target
%   probability
% 
%   C - Bayesian coverage goal
%   
%   ALPHA, BETA - parameters of Beta prior
%
%   TRACES_RES - the results of model-checking traces (0: fail, 1: success)
%
%
% $Revision: 1.2.1 $
% $Date: 2012/12/19 14:07:32 $
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n = 0;   % initial total number of samples
x = 0;   % initial number of success samples
r = 0;   % initial posterior probability of estimation interval
cnt = 1; % initial index of traces_res

while r < c % until r satisfies user-given coverage goal
  % add 1 to total number of samples
  n = n + 1;
  
  % if cnt's result of model-checking traces is success,
  % add 1 to number of success samples
  if traces_res(cnt) == 1 
      x = x + 1;
  end
  cnt = cnt + 1;
  
  % compute the Beta prior probability with alpha and beta
  p = (x + alpha)/(n + alpha + beta);
  
  % compute the estimation interval which may contain target probability
  t0 = p - delta;
  t1 = p + delta;
  if t1 > 1
      t0 = 1 - 2*delta;
      t1 = 1;
  elseif t0 < 0
      t0 = 0;
      t1 = 2*delta;
  end
  
  % compute posterior probability r that target probability
  % is included within interval (t0,t1)
  r = posterior_prob(t0,t1,n,x,alpha,beta);
  
  % log for each step
  str=sprintf('t0: %f t1: %f p: %f n: %d x: %d\n',t0,t1,p,n,x);
  disp(str);
end

end

function result = posterior_prob(t0,t1,n,x,alpha,beta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% posterior_prob - Function of calculating the posterior probability of
% Bernoulli samples and Beta prior of parameters alpha and beta
%
%  POSTERIOR_PROB(T0,T1,N,X,ALPHA,BETA)
%
%  T0 - left-open value of estimated interval
%  T1 - right-open value of estimated interval
%  N - total number of samples
%  X - number of success Bernoulli samples
%  ALPHA, BETA - parameters of Beta prior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
result = betacdf(t1,x+alpha,n-x+beta)-betacdf(t0,x+alpha,n-x+beta);
end
