diff --git a/examples/13-genint/main.cpp b/examples/13-genint/main.cpp index 0f8d039c..54149240 100644 --- a/examples/13-genint/main.cpp +++ b/examples/13-genint/main.cpp @@ -10,53 +10,14 @@ int main() { // Model parameters - double r = 5.0; // Contact rate - double popsize = 1e4; - double p_contact = r/popsize; // Probability of contact - double p_i = .1; // Probability of infection - double p_r = 1.0/15.0; // Probability of recovery - - size_t nsims = 1e6; // Number of simulations - size_t nsteps = 100; // Number of days - - // Generation intervals: time to first transmission - std::vector< double > genints; - genints.reserve(nsims); - - // Random number generator - std::mt19937 rng; - std::uniform_real_distribution unif(0.0, 1.0); - std::binomial_distribution binom(popsize, p_contact); - rng.seed(0); - - for (size_t i = 0; i < nsims; i++) - { - - for (size_t j = 0; j < nsteps; j++) { - - // Recovers - if (unif(rng) < p_r) - break; - - // How many contacts - int ncontacts = binom(rng); - - // Infects at least 1 - double p = 1.0 - pow(1.0 - p_i, ncontacts); - if (unif(rng) < p) { - genints.push_back(j); - break; - } - } - } - + double r = 10.0; // Contact rate + double popsize = 1e5; + double p_i = .1; // Probability of infection + double p_r = 1.0/7.0; // Probability of recovery + double incubation = 2.0; // Incubation period - // Computing the mean - double mean = 0.0; - for (auto & genint : genints) - mean += genint; - - mean /= genints.size(); + size_t nsteps = 200; // Number of days + size_t max_days_for_calc = 20; epiworld::epimodels::ModelSIRCONN<> model( "avirus", @@ -70,23 +31,30 @@ int main() model.run(nsteps, 554); model.print(); - std::vector agent, virus, time, gentime; - model.get_db().generation_time(agent, virus, time, gentime); - + std::vector agent, virus, time, gentime_sir_obs; + model.get_db().generation_time(agent, virus, time, gentime_sir_obs); + auto gentime_sir = model.generation_time_expected(); + double gentime_sir_mean = std::accumulate( + gentime_sir.begin(), gentime_sir.begin() + max_days_for_calc, 0.0 + ) / max_days_for_calc; + // Computing the mean of gentime conditioning on gentime >= 0 - double mean2 = 0.0; + double gentime_sir_obs_mean = 0.0; size_t n = 0; - for (int i = 0; i < gentime.size()*3/4; i++) + for (int i = 0; i < gentime_sir_obs.size()*3/4; i++) { - if (gentime[i] >= 0) + if (time[i] >= max_days_for_calc) + continue; + + if (gentime_sir_obs[i] >= 0) { - mean2 += gentime[i]; + gentime_sir_obs_mean += gentime_sir_obs[i]; n++; } } - mean2 /= n; + gentime_sir_obs_mean /= n; // Repeating for ModelSEIRCONN epiworld::epimodels::ModelSEIRCONN<> model2( @@ -95,38 +63,43 @@ int main() 50.0/popsize, r, p_i, - 0.00000000009, + incubation, p_r ); model2.run(nsteps, 554); model2.print(); - std::vector agent2, virus2, time2, gentime2; - model2.get_db().generation_time(agent2, virus2, time2, gentime2); + std::vector agent2, virus2, time2, gentime_seir_obs; + model2.get_db().generation_time(agent2, virus2, time2, gentime_seir_obs); + auto gentime_seirconn = model2.generation_time_expected(); + double gentime_seirconn_mean = std::accumulate( + gentime_seirconn.begin(), gentime_seirconn.begin() + max_days_for_calc, 0.0 + ) / max_days_for_calc; // Computing the mean of gentime conditioning on gentime >= 0 - double mean_seirconn = 0.0; + double gentime_seir_obs_mean = 0.0; n = 0; - for (int i = 0; i < gentime2.size()*3/4; i++) + for (int i = 0; i < gentime_seir_obs.size()*3/4; i++) { - if (gentime2[i] >= 0) + if (time[i] >= max_days_for_calc) + continue; + if (gentime_seir_obs[i] >= 0) { - mean_seirconn += gentime2[i]; + gentime_seir_obs_mean += gentime_seir_obs[i]; n++; } } - mean_seirconn /= n; + gentime_seir_obs_mean /= n; - // Printing both mean and mean2 - std::cout << "Mean (naive) : " << mean << std::endl; - std::cout << "Mean2 (from sim) : " << mean2 << std::endl; - std::cout << "Mean (SEIRCONN) : " << mean_seirconn << std::endl; - double theory_mean = epiworld::gen_int_mean(popsize, p_contact, p_i, p_r); - std::cout << "Mean (theory) +1: " << theory_mean + 1.0 << std::endl; + std::cout << "SIR Gen. Int. (obs) : " << gentime_sir_obs_mean << std::endl; + std::cout << "SIR Gen. Int. (expected) : " << gentime_sir_mean << std::endl; + + std::cout << "SEIR Gen. Int. (obs) : " << gentime_seir_obs_mean << std::endl; + std::cout << "SEIR Gen. Int. (expected) : " << gentime_seirconn_mean << std::endl; return 0; diff --git a/include/epiworld/math/distributions.hpp b/include/epiworld/math/distributions.hpp index 0460e70a..b32a012e 100644 --- a/include/epiworld/math/distributions.hpp +++ b/include/epiworld/math/distributions.hpp @@ -76,7 +76,8 @@ inline double dgenint( int max_days = 200 ) { - g += 1; + if ((g < 1) || (g > max_days)) + return 0.0; if (p_0_approx < 0.0) { @@ -155,7 +156,7 @@ inline double gen_int_mean( mean += static_cast(i) * dgenint( - i, S, p_c, p_i, p_r, p_0_approx, normalizing, max_n + i, S, p_c, p_i, p_r, p_0_approx, normalizing, max_n, max_days ); } diff --git a/include/epiworld/models/seirconnected.hpp b/include/epiworld/models/seirconnected.hpp index 999dfc43..296f2f3c 100644 --- a/include/epiworld/models/seirconnected.hpp +++ b/include/epiworld/models/seirconnected.hpp @@ -431,16 +431,16 @@ inline std::vector< double > ModelSEIRCONN::generation_time_expected( // spend at least one day in the exposed state, and 1 day in the // infectious state before starting transmitting. std::vector< double > gen_times( - this_const->get_ndays(), 2.0 + days_exposed + this_const->get_ndays(), 1.0 + days_exposed ); double p_c = this_const->par("Contact rate")/this_const->size(); - double p_i = this_const->par("Transmission rate"); - double p_r = this_const->par("Recovery rate"); + double p_i = this_const->par("Prob. Transmission"); + double p_r = this_const->par("Prob. Recovery"); for (size_t i = 0u; i < this_const->get_ndays(); ++i) { - gen_times[i] = gen_int_mean( + gen_times[i] += gen_int_mean( S[i], p_c, p_i, diff --git a/tests/10-generation-interval.cpp b/tests/10-generation-interval.cpp index e9178c9c..1162ecc7 100644 --- a/tests/10-generation-interval.cpp +++ b/tests/10-generation-interval.cpp @@ -14,8 +14,8 @@ EPIWORLD_TEST_CASE("Generation interval", "[gen-int]") static const int S = 10000; static const size_t max_days = 200; static const size_t max_contacts = S/2; - static const double p_r = 1.0/15.0; - static const double p_i = 0.1; + static const double p_r = 1.0/7.0; + static const double p_i = 0.05; // Derived parameters static const double p_c = contact_rate/static_cast(S); @@ -63,7 +63,7 @@ EPIWORLD_TEST_CASE("Generation interval", "[gen-int]") if (model.runif() < p_at_least_one) { p_0 -= 1.0; - sim_days.push_back(j); + sim_days.push_back(j + 1); distribution[j] += 1.0; break; } @@ -95,7 +95,7 @@ EPIWORLD_TEST_CASE("Generation interval", "[gen-int]") for (size_t i = 0u; i < (max_days * 5); ++i) { distribution_expected[i] = epiworld::dgenint( - i, S, p_c, p_i, p_r, p_0_approx_analy_global, + i + 1.0, S, p_c, p_i, p_r, p_0_approx_analy_global, normalizing_constant, max_contacts, max_days ); }