Add files via upload

Examples/Example1_BellmanHarrisExpectedPopulation.m

@@ -0,0 +1,73 @@
+%% Example 1: Calculating the expected future population in Belman-Harris Branching Process (BHBP)
+% The most used distribution for the life length in BHBP is exponential distribution, because it is the only case 
+% in which we know the theoretical expectation for the population count. In all other cases we know that it satisfies 
+% a certain renewal equation (Z(t) = S + \int_0^t [m*Z(t-y)] dG(y)) but we do not know its theoretical solution.
+% This example considers the exponential life length and compares the numerical and theoretical solution. 
+% Then it compares the results for exponential and truncated normal life length.
+%
+%       Copyright (C) 2018  Plamen Ivaylov Trayanov
+%
+%       This program is free software: you can redistribute it and/or modify
+%       it under the terms of the GNU General Public License as published by
+%       the Free Software Foundation, either version 3 of the License, or
+%       (at your option) any later version.
+% 
+%       This program is distributed in the hope that it will be useful,
+%       but WITHOUT ANY WARRANTY; without even the implied warranty of
+%       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%       GNU General Public License for more details.
+% 
+%       You should have received a copy of the GNU General Public License
+%       along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%
+addpath('../')      % adds the function numerical_solution to the path
+T=120;      % defines the interval [0, T] in which we want to find numerical solution, i.e. the expected population count
+h=0.1;      % defines the grid size of the numerical method
+
+mean_life_length=10;    % define the average life length to be 10
+S=1-expcdf(0:h:T, mean_life_length)';   % defines the survivability function of the cells
+H=[3/4,0,1/4]';     % defines the p.g.f of the offspring: P(0 offpsring) = 3/4, P(2 offspring) = 1/4, or in other words
+% the p.g.f. is 3/4 + 0*s^1 + 1/4*s^2. 
+
+% defines the functions used in the integral equation
+G=1-S;      % the life length distribution
+m=(0:(length(H)-1))*H;      % m=0.5 so the process is subcritical.
+f=@(y, Z)(m*Z);     % the integrand in the equation
+z=S;
+
+Z_numerical_exp=numerical_solution(z, f, G, T, h);    % finds the solution, i.e. the expected population count in a Bellman-Harris branching process
+Z_true=exp((m-1)*(0:h:T)/mean_life_length);     % the theoretical solution
+
+% calculates the solution for normally distributed life length (truncated normal as the life length cannot be negative)
+% here we do not have theoretical distribution, so numerical or simulation method is the only way to solve this
+S=(1-normcdf(0:h:T, mean_life_length, 2.5)')./(1-normcdf(0, mean_life_length, 2.5));  % defines the survivability function of the cells
+G=1-S;
+z=S;
+Z_numerical_normal=numerical_solution(z, f, G, T, h);
+
+%% make plots of the results
+% plot the numerical vs theoretical solution when the life length is exponentially distributed
+line_wd=2.5;
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:T, Z_numerical_exp, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_true, '--', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Expected population count');
+legend('Numerical solution', 'Theoretical Solution')
+title('Bellman-Harris BP with exponential life length Exp(10)')
+print('./figures/Example1_fig1', '-dpng', '-r0')
+
+% plot the expected population count fo exponentially distributed vs normally distributed life length
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:T, Z_numerical_exp, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_numerical_normal, '--', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Expected population count');
+legend('Exponential life length', 'Normal life length')
+title('Bellman-Harris BP with mean life length of 10')
+print('./figures/Example1_fig2', '-dpng', '-r0')
+
+

Examples/Example2_SuccessfulMutants.m

@@ -0,0 +1,113 @@
+%% Example 2: Calculating the distribution of successful mutant arrival in two-type Belman-Harris Branching Process (BHBP)
+% This example considers the model in "Branching processes in continuous time as models of mutations: Computational 
+% approaches and algorithms by Slavtchova-Bojkova, M., Trayanov, P., Dimitrov, S. (2017). Computational Statistics and 
+% Data Analysis, http://dx.doi.org/10.1016/j.csda.2016.12.013.
+%
+% Short description of the article model:
+% type 0 particles - supercritical, do not mutate to other types once they occur
+% type 1 particles - subcritical, with probability of mutation to type 0
+% The population starts from 1 particle of type 1.
+% "Successful mutant" - a particle of type 0 that starts a type 0 branching, that never extincts
+% "Unsuccessful mutant" - a particle of type 0 that starts a type 0 branching, that extincts
+% The object of interest in the article is to calculate the distribution T of the first arrival of a "successful mutant".
+%
+% The distribution of the arrival time satisfies the integral equation:
+% P(T>t) = 1-G(t) + \int_0^t f_1(u*q_0+(1-u)*P(T>t-y)) dG(y)
+% The example shows the results for different distributions of life lengths.
+%
+%       Copyright (C) 2018  Plamen Ivaylov Trayanov
+%
+%       This program is free software: you can redistribute it and/or modify
+%       it under the terms of the GNU General Public License as published by
+%       the Free Software Foundation, either version 3 of the License, or
+%       (at your option) any later version.
+% 
+%       This program is distributed in the hope that it will be useful,
+%       but WITHOUT ANY WARRANTY; without even the implied warranty of
+%       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%       GNU General Public License for more details.
+% 
+%       You should have received a copy of the GNU General Public License
+%       along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%
+addpath('../')      % adds the function numerical_solution to the path
+T=100;      % defines the interval [0, T] in which we want to find numerical solution, i.e. the expected population count
+h=0.1;      % defines the grid size of the numerical method
+
+mean_life_length=10;    % we want to have average life length of 10
+S=1-expcdf(0:h:T, mean_life_length)';   % defines the survivability function of the cells
+H=[1-0.375, 0, 0.375]';      % defines the p.g.f of the offspring: P(0 offpsring) = 3/4, P(2 offspring) = 1/4, or in other words
+% the p.g.f. is 3/4 + 0*s^1 + 1/4*s^2. 
+
+q_0=0.3;    % extinction probability for type 0
+u=0.2;      % mutation probability from type 1 to type 0
+
+% defines the functions used in the integral equation
+G=1-S;      % the life length distribution
+z=S;
+f_1=@(x)(arrayfun(@(x)(sum(H'.*x.^(0:size(H,1)-1))), x));   % the offspring p.g.f. for type 1 cells
+f=@(y, Z)(f_1(u*q_0+(1-u)*Z));      % defines the integrand in the integral equation
+
+Z_exp=numerical_solution(z, f, G, T, h);    % finds the solution, i.e. P(T>t), for exponential distribution G(t)
+
+% calculates the solution for normally distributed life length
+S=(1-normcdf(0:h:T, mean_life_length, 2.5)')./(1-normcdf(0, mean_life_length, 2.5));  % defines the survivability function of the cells
+Z_normal=numerical_solution(S, f, 1-S, T, h);   % finds the solution for truncated normal life length G(t)
+
+% calculates the solution for uniformly distributed life length
+S=1-unifcdf(0:h:T, 0, mean_life_length*2)';  % defines the survivability function of the cells
+Z_unif=numerical_solution(S, f, 1-S, T, h);
+
+% calculates the solution for beta distributed life length
+S=[1-betacdf(0:h/(2*mean_life_length):1, 2, 2)'; zeros(T/h-length(0:h/(2*mean_life_length):1)+1, 1)];
+Z_beta=numerical_solution(S, f, 1-S, T, h);
+
+% calculates the solution for chi^2 distributed life length
+S=1-chi2cdf(0:h:T, mean_life_length)';
+Z_chi2=numerical_solution(S, f, 1-S, T, h);
+
+%% make plots of the results
+line_wd=2.5;
+% plot the different life length distributions used in the example
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:T, exppdf(0:h:T, 20)*h, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, normpdf(0:h:T, 10, 2.5)*h, '-', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, unifpdf(0:h:T, 0, 20)*h, '-', 'Color', [0, 1, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:20, betapdf(0:(h/20):1, 2, 2)*h/20, '-', 'Color', [0, 0, 1, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, chi2pdf(0:h:T, 10)*h, '--', 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Probability');
+legend('Exp(10) life length', 'N(10, 2.5) life length', 'U(0, 20) life length', 'Beta(2,2) life length, scaled in [0, 20]', 'Chi^2(10) life length')
+title('Life length distributions used in the example')
+print('./figures/Example2_fig1', '-dpng', '-r0')
+
+% plot the probability that successful mutant has not arrived yet for different life length distributions
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:T, Z_exp, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_normal, '-', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_unif, '-', 'Color', [0, 1, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_beta, '-', 'Color', [0, 0, 1, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_chi2, '--', 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Probability P(T>t)');
+legend('Exp(10) life length', 'N(10, 2.5) life length', 'U(0, 20) life length', 'Beta(2,2) life length, scaled in [0, 20]', 'Chi^2(10) life length')
+title('Successful mutant arrival in two-type Bellman-Harris BP')
+print('./figures/Example2_fig2', '-dpng', '-r0')
+
+% plot the distribution density function of T (time of successful mutant arrival)
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(h:h:T, -diff(Z_exp)*h, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(h:h:T, -diff(Z_normal)*h, '-', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(h:h:T, -diff(Z_unif)*h, '-', 'Color', [0, 1, 0, 0.5], 'LineWidth', line_wd)
+plot(h:h:T, -diff(Z_beta)*h, '-', 'Color', [0, 0, 1, 0.5], 'LineWidth', line_wd)
+plot(h:h:T, -diff(Z_chi2)*h, '--', 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Probability Density Function of T');
+legend('Exponential', 'Normal', 'Uniform', 'Beta', 'ChiSquared')
+title('Changing the life length distribution')
+print('./figures/Example2_fig3', '-dpng', '-r0')
+

Examples/Example3_SuccessfulMutants.m

@@ -0,0 +1,78 @@
+%% Example 3: Calculating the distribution of successful mutant arrival in two-type Belman-Harris Branching Process (BHBP)
+% This example is continuation of the previous Example2. Here we discuss the joint probability that "successful mutant" 
+% has not been born until time t and we do not have cells of type 1 (with subcritical reproduction rate, m_1 < 1) at 
+% time t. This joint probability satisfies an integral equation that we will solve for different life length 
+% distributions.
+% P(T>t, Z^1(t)=0)=\int_0^t f_1(uq_0+(1-u)P(T>t-y, Z^1(t-y)=0)) dG_1(y).
+%
+%
+%       Copyright (C) 2018  Plamen Ivaylov Trayanov
+%
+%       This program is free software: you can redistribute it and/or modify
+%       it under the terms of the GNU General Public License as published by
+%       the Free Software Foundation, either version 3 of the License, or
+%       (at your option) any later version.
+% 
+%       This program is distributed in the hope that it will be useful,
+%       but WITHOUT ANY WARRANTY; without even the implied warranty of
+%       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%       GNU General Public License for more details.
+% 
+%       You should have received a copy of the GNU General Public License
+%       along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%
+addpath('../')      % adds the function numerical_solution to the path
+T=100;      % defines the interval [0, T] in which we want to find numerical solution, i.e. the expected population count
+h=0.1;      % defines the grid size of the numerical method
+
+mean_life_length=10;    % we want to have average life length of 10
+S=1-expcdf(0:h:T, mean_life_length)';   % defines the survivability function of the cells
+H=[1-0.375, 0, 0.375]';      % defines the p.g.f of the offspring: P(0 offpsring) = 1-0.375, P(2 offspring) = 0.375, 
+% or in other words the p.g.f. is 1-0.375 + 0*s^1 + 0.375*s^2. 
+
+q_0=0.3;    % extinction probability
+u=0.2;      % mutation probability
+
+% defines the functions used in the integral equation
+G=1-S;      % the life length distribution
+z=zeros(size(S));
+f_1=@(x)(arrayfun(@(x)(sum(H'.*x.^(0:size(H,1)-1))), x));   % the offspring p.g.f. for type 1 cells
+f=@(y, Z)(f_1(u*q_0+(1-u)*Z));      % the equation as in article http://dx.doi.org/10.1016/j.csda.2016.12.013
+
+Z_exp=numerical_solution(z, f, G, T, h);    % finds the solution
+
+% calculates the solution for normally distributed life length
+S=(1-normcdf(0:h:T, mean_life_length, 2.5)')./(1-normcdf(0, mean_life_length, 2.5));  % defines the survivability function of the cells
+Z_normal=numerical_solution(z, f, 1-S, T, h);
+
+% calculates the solution for uniformly distributed life length
+S=1-unifcdf(0:h:T, 0, mean_life_length*2)';  % defines the survivability function of the cells
+Z_unif=numerical_solution(z, f, 1-S, T, h);
+
+% calculates the solution for beta distributed life length
+S=[1-betacdf(0:h/(2*mean_life_length):1, 2, 2)'; zeros(T/h-length(0:h/(2*mean_life_length):1)+1, 1)];  % defines the survivability function of the cells
+Z_beta=numerical_solution(z, f, 1-S, T, h);
+
+% calculates the solution for chi^2 distributed life length
+S=1-chi2cdf(0:h:T, mean_life_length)';  % defines the survivability function of the cells
+Z_chi2=numerical_solution(z, f, 1-S, T, h);
+
+
+%% make plots of the results
+line_wd=2.5;
+% plot the probability that successful mutant has not arrived yet
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:T, Z_exp, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_normal, '-', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_unif, '-', 'Color', [0, 1, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_beta, '-', 'Color', [0, 0, 1, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, Z_chi2, '--', 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Probability P(T>t, Z^1(t)=0)');
+legend('Exp(10) life length', 'N(10, 2.5) life length', 'U(0, 20) life length', 'Beta(2,2) life length, scaled in [0, 20]', ...
+    'Chi^2(10) life length', 'Location', 'SE')
+title('Successful mutant arrival in two-type Bellman-Harris BP')
+print('./figures/Example3_fig1', '-dpng', '-r0')
+

Examples/Example4_SolvingRenewalEquations.m

@@ -0,0 +1,91 @@
+%% Example 4: Solving Renewal Equations and estimating the renewal functions
+% A renewal equation is integral equation of the type Z(t) = z(t) + \int_0^t Z(t-y) dG(y).
+% In Renewal Theory its solution is represented by infinite sum of convolutions of increasing order:
+% Z(t) = U*z(t), where U is the renewal function G^(0*)+ G^(1*) + ... + G^(n*) + ...
+% This numerical method however does not calculate these convolutions numerically and is much faster.
+% Also, the renewal function satisfies the equation U(t) = I(t) + \int_0^t U(t-y) dG(t) and we can find its solution 
+% without calculating the convolutions. U(t)=N(t)+1, the number of renewals, counting the initial moment 0 as a renewal
+% time.
+%
+%
+%       Copyright (C) 2018  Plamen Ivaylov Trayanov
+%
+%       This program is free software: you can redistribute it and/or modify
+%       it under the terms of the GNU General Public License as published by
+%       the Free Software Foundation, either version 3 of the License, or
+%       (at your option) any later version.
+% 
+%       This program is distributed in the hope that it will be useful,
+%       but WITHOUT ANY WARRANTY; without even the implied warranty of
+%       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%       GNU General Public License for more details.
+% 
+%       You should have received a copy of the GNU General Public License
+%       along with this program. If not, see <http://www.gnu.org/licenses/>.
+%
+%
+addpath('../')      % adds the function numerical_solution to the path
+T=100;          % defines the interval [0, T] in which we want to find numerical solution, i.e. the expected population count
+h=0.01;         % defines the grid size of the numerical method
+
+%% define a Poisson renewal process with averate rate 1/10
+% the expected number of renewals (without counting the initial moment 0 as a renewal time) satisfies the equation
+% N(t) = G(t) + \int_0^t N(t-y) dG(y).
+% For Poisson process we know the solution is actually Z(t)=lambda*t, where lambda is the renewal rate (1/<the average 
+% renewal interval length>). We can compare the theoretical and numerical solutions.
+% The renewal function U(t) satisfies U(t) = I(t) + U * G
+
+mean_life_length=10;    % we want to have an average renewal time of 10
+G=expcdf(0:h:T, mean_life_length)';   % defines the distributions of the time between renewals to be exponential
+z=G;
+f=@(y, Z)(Z);   % defines the equation as renewal-type equation
+
+% the expected number of renewals satisfies this equation
+N_poisson_numerical=numerical_solution(z, f, G, T, h);    % estimates the solution, i.e. the expected number of renewals
+N_poisson_theoretical=(0:h:T)/mean_life_length;
+
+% the renewal function for Poisson process:
+z=ones(T/h+1,1);
+U_poisson_numerical=numerical_solution(z, f, G, T, h);  % estimates the renewal function as a solution to renewal equation
+U_poisson_theoretical=1+(0:h:T)/mean_life_length;       % the theoretical solution is known in this example
+
+% make plots of the results
+line_wd=2.5;
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:T, N_poisson_numerical, 'Color', [1, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, N_poisson_theoretical, '--', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, U_poisson_numerical, 'Color', [0, 0, 0.7, 0.5], 'LineWidth', line_wd)
+plot(0:h:T, U_poisson_theoretical, '--', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Expected number of renewals');
+legend({'Numerical solution for $E(N(t))$', 'Theoretical Solution: $E(N(t))=\lambda t$', 'Numerical solution for $U(t)$', ...
+    'Theoretical Solution: $U(t)=1+\lambda t$'}, 'interpreter', 'latex')
+title('Poisson process. $N(t) = G(t) + N\ast G (t), U(t) = I(t) + U\ast G (t)$','interpreter', 'latex')
+print('./figures/Example4_fig1', '-dpng', '-r0')
+
+%% find the renewal function for uniform distribution G in [0, 1]
+% in this case we know the theoretical solution for t in the interval [0, 1]: U(t) = e^t
+T=1;
+h=0.001;
+G=unifcdf(0:h:T, 0, 1)';    % defines the distributions of the time between renewals to be exponential
+f=@(y, Z)(Z);               % defines the equation as renewal type equation
+
+z=ones(T/h+1,1);
+U_uniform_numerical=numerical_solution(z, f, G, T, h);
+U_uniform_theoretical=exp(0:h:T);
+
+line_wd=2.5;
+figure('visible','on', 'Units','pixels','OuterPosition',[0 0 1280 1024]);
+set(gca,'FontSize',16)
+hold on
+plot(0:h:1, U_uniform_numerical, 'Color', [0, 0, 0.7, 0.5], 'LineWidth', line_wd)
+plot(0:h:1, U_uniform_theoretical, '--', 'Color', [0, 0, 0, 0.5], 'LineWidth', line_wd)
+xlabel('Time'); ylabel('Renewal Function');
+legend({'Numerical solution for $U(t)$', 'Theoretical solution: $U(t)=e^t$ for $t \in [0,1]$'}, ...
+    'Location', 'NW', 'interpreter', 'latex')
+title('G is defined to be uniform in [0, 1]')
+print('./figures/Example4_fig2', '-dpng', '-r0')
+
+
+