Updating examples

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<double> unif(0.0, 1.0);
-    std::binomial_distribution<int> 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<int> agent, virus, time, gentime;
-    model.get_db().generation_time(agent, virus, time, gentime);
-
+    std::vector<int> 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<int> agent2, virus2, time2, gentime2;
-    model2.get_db().generation_time(agent2, virus2, time2, gentime2);
+    std::vector<int> 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;

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<double>(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
                 );
 
     }

include/epiworld/models/seirconnected.hpp

@@ -431,16 +431,16 @@ inline std::vector< double > ModelSEIRCONN<TSeq>::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,

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<double>(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
         );
     }