% % 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%