% % CMC_main.m % % Use crude Monte Carlo (indpendent sample paths) 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: CMC_main(n,tend,theta,epsilon) % % n = number of times to run CMC_path to get stats. % % tend = terminal time % % theta = value of parameter at which we are computing the % derivative. % % epsilon = epsilon used in centered finite difference % % 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: CMC_main(5000,3,0.5,0.5*1/20) % function [m_C,m_d,v_d] = CMC_main(n,tend,theta,epsilon) %initial condition initial = [300; 0; 0]; %"data_plus" will store the output with theta + epsilon/2 at time tend. %Each column will represent a realization. data_plus = zeros([3,n]); %"data_minus" will store the output with theta - epsilon/2 at time tend. %Each column will represent a realization. data_minus = zeros([3,n]); %"derivative" will give the estimated derivative achieved by taking a %difference and dividing by epsilon. derivative = zeros([n,1]); for j = 1:n [data_plus(:,j)] = CMC_path(tend,initial,theta + epsilon/2); [data_minus(:,j)] = CMC_path(tend,initial,theta - epsilon/2); derivative(j) = (data_plus(3,j) - data_minus(3,j))/epsilon; end m_C = mean(data_plus(3,:) + data_minus(3,:))/2; m_d = mean(derivative); v_d = var(derivative)/n; end %The program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%