%
%  Likelihood_main.m
%
%           Use Likelihood ratio method to compute the sensitivity of the
%           expected number of C molecules with respect to theta for the 
%           model
%
%                   A -> B -> C  
%
%               with rates k1 = theta, k2 = 1.  
%
%  usage:	Likelihood_main(n,tend,theta)
%
%               n = number of times to run Likelihood_path to get stats.
%
%               tend = terminal time
%
%               theta = value of parameter at which we are computing the
%                       derivative.
%
%  output:  m_C, m_d, v_d
%
%               m_C = estimated value of E C(tend).
%               m_d = estimated value of the derivative of E C(tend)
%               v_d = variance of estimate for m_d.
%
%  example:	Likelihood_main(5000,3,0.5)
%

function [m_C,m_d,v_d] = Likelihood_main(n,tend,theta)

%initial condition
initial = [300; 
            0; 
            0];

%"data" will store the output at time tend.  Each column will represent a
%realization.
data = zeros([3,n]);

%"M" stores the weights used to compute the derivatives.
M = zeros([n,1]);

%"derivative" will give the estimated derivative achieved by multiplying
%data(3,j) by M(j).
derivative = zeros([n,1]);

%Call Likelihood_path n times to get the statistics.
for j = 1:n
  
    %data(:,j) gives the jth output.  M gives the appropriate weighting
    %function.
    [data(:,j),M(j)] = Likelihood_path(tend,initial,theta);
  
    %the estimated derivative.
    derivative(j) = data(3,j)*M(j);
      
end

m_C = mean(data(3,:));
m_d = mean(derivative);
v_d = var(derivative)/n;


end %The program

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%