fix error in filtering of SEIR and Two-Species models

Author: kingaa

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: phylopomp
 Type: Package
 Title: Phylodynamic Inference for POMP Models
-Version: 0.14.7.0
-Date: 2024-11-06
+Version: 0.14.8.0
+Date: 2024-11-15
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="[email protected]",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Qianying"),family="Lin",role=c("aut"),comment=c(ORCID="0000-0001-8620-9910"))
 	  )

inst/NEWS

@@ -1,14 +1,15 @@
 _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_h_y_l_o_p_o_m_p'
 
-_C_h_a_n_g_e_s _i_n '_p_h_y_l_o_p_o_m_p' _v_e_r_s_i_o_n _0._1_6:
+_C_h_a_n_g_e_s _i_n '_p_h_y_l_o_p_o_m_p' _v_e_r_s_i_o_n _0._1_4:
 
         • The filtering codes for the two-species model are now
           working.
 
         • A new sample/death event-type has been introduced into the
           ‘genealogy_t’ type.
 
-_C_h_a_n_g_e_s _i_n '_p_h_y_l_o_p_o_m_p' _v_e_r_s_i_o_n _0._1_4:
+        • Corrections of filtering equations for the SEIR and
+          Two-species models.
 
         • Yet another refactoring.  The SEIRS, LBDP, and SIRS models
           have been refactored to streamline the filtering calculations

inst/NEWS.Rd

@@ -1,13 +1,10 @@
 \name{NEWS}
 \title{News for package `phylopomp'}
-\section{Changes in \pkg{phylopomp} version 0.16}{
+\section{Changes in \pkg{phylopomp} version 0.14}{
   \itemize{
     \item The filtering codes for the two-species model are now working.
     \item A new sample/death event-type has been introduced into the \code{genealogy_t} type.
-  }
-}
-\section{Changes in \pkg{phylopomp} version 0.14}{
-  \itemize{
+    \item Corrections of filtering equations for the SEIR and Two-species models.
     \item Yet another refactoring.
     The SEIRS, LBDP, and SIRS models have been refactored to streamline the filtering calculations in a generalizable fashion.
   }

src/seirs_pomp.c

@@ -288,10 +288,10 @@ void seirs_gill
         break;
       case 2:                   // transmission, s=(1,0)
         assert(S>=1);
-        S -= 1; E += 1;
-        ll += log(1-ellI/I)-log(E);
         change_color(color,nsample,random_choice(ellI),1,0);
         ellE += 1; ellI -= 1;
+        S -= 1; E += 1;
+        ll += log(1-ellI/I)-log(E);
         assert(!ISNAN(ll));
         break;
       case 3:                   // progression, s=(0,0)
@@ -302,10 +302,10 @@ void seirs_gill
         break;
       case 4:                   // progression, s=(0,1)
         assert(E>=1);
-        E -= 1; I += 1;
-        ll -= log(I);
         change_color(color,nsample,random_choice(ellE),0,1);
         ellE -= 1; ellI += 1;
+        E -= 1; I += 1;
+        ll -= log(I);
         assert(!ISNAN(ll));
         break;
       case 5:                   // recovery

src/twospecies_pomp.c

@@ -114,13 +114,17 @@ static double event_rates
   *logpi = 0; logpi++;
   assert(R_FINITE(event_rate));
   // 2: Trans_21, s = (0,0),(1,0)
+  // s = (0,0): pi = (I1-ell1)/I1, m = 1,
+  // s = (1,0): pi = 1/(2*I1), m = ell1
+  // total pi = (I1-0.5*ell1)/I1
   assert(S2>=0 && I1>=ell1 && ell1>=0);
   alpha = (N1 > 0) ? Beta21*S2*I1/N1 : 0;
   pi = (I1 > 0) ? (I1-0.5*ell1)/I1 : 0;
   event_rate += (*rate = alpha*pi); rate++;
   *logpi = log(pi); logpi++;
   assert(R_FINITE(event_rate));
   // 3: Trans_21, s = (0,1)
+  // s = (0,1): pi = 1/(2*I1), m = ell1
   pi = 1-pi;
   event_rate += (*rate = alpha*pi); rate++;
   *logpi = log(pi)-log(ell1); logpi++;
@@ -531,9 +535,9 @@ void twospecies_gill
       case 3:                   // 3: Trans_21, s = (0,1)
         assert(S2>=1 && I1>=0);
         S2 -= 1; I2 += 1;
-        ll += log(1-ell1/I1)-log(I2);
         change_color(color,nsample,random_choice(ell1),host1,host2);
         ell1 -= 1; ell2 += 1;
+        ll += log(1-ell1/I1)-log(I2);
         assert(check_color(color,nsample,ell1,ell2));
         assert(!ISNAN(ll));
         break;
@@ -546,9 +550,9 @@ void twospecies_gill
       case 5:                   // 5: Trans_12, s = (1,0)
         assert(S1>=1 && I2>=0);
         S1 -= 1; I1 += 1;
-        ll += log(1-ell2/I2)-log(I1);
         change_color(color,nsample,random_choice(ell2),host2,host1);
         ell2 -= 1; ell1 += 1;
+        ll += log(1-ell2/I2)-log(I1);
         assert(check_color(color,nsample,ell1,ell2));
         assert(!ISNAN(ll));
         break;

tests/curtail.Rout.save

@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

tests/diagram.Rout.save

@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

tests/gendat.Rout.save

@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

tests/lbdp1.Rout.save

@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu
 

tests/lbdp2.Rout.save

@@ -1,5 +1,5 @@
 
-R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+R version 4.4.2 (2024-10-31) -- "Pile of Leaves"
 Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu