two bugfixes

.Rproj.user/BAB4038F/pcs/workbench-pane.pper

@@ -1,5 +1,5 @@
 {
     "TabSet1": 1,
-    "TabSet2": 5,
+    "TabSet2": 0,
     "TabZoom": {}
 }

R/aoy.R

@@ -46,7 +46,7 @@ dAoY = function(alpha, a, FoIpar, tau=0, hhat=1, r=1/200){
   # The function call
   dAoYcompute = function(alpha, a, FoIpar, tau, hhat, r){
     moi = meanMoI(a,FoIpar,tau,hhat,r)
-    moi*dAoI(alpha,a,FoIpar,tau,hhat,r)*exp(-moi*pAoI(alpha,a,FoIpar,tau,hhat,r))/truePRa(a,FoIpar,hhat,tau,r)
+    moi*dAoI(alpha,a,FoIpar,tau,hhat,r)*exp(-moi*pAoI(alpha,a,FoIpar,tau,hhat,r))/truePRa(a,FoIpar,tau,hhat,r)
   }
   # Use sapply to call dAoYcompute multiple times
   if(length(alpha)==1) return(dAoYcompute(alpha, a,FoIpar,tau,hhat,r))
@@ -86,7 +86,7 @@ momentAoY = function(a, FoIpar, tau=0, hhat=1, r=1/200, n=1){
     ff = function(alpha, a, FoIpar, tau, hhat, r, n){
       alpha^n*dAoY(alpha, a, FoIpar, tau, hhat, r)
     }
-    stats::integrate(ff, 0, a,a=a,FoIpar=FoIpar, tau=tau, hhat=hhat,r=r,n=n)$value
+    stats::integrate(ff, 0, a,a=a, FoIpar=FoIpar, tau=tau, hhat=hhat,r=r,n=n)$value
   }
   if(length(a)==1){return(ffAoYda(a, FoIpar, tau, hhat, r, n))} else{
     sapply(a, ffAoYda, FoIpar=FoIpar, tau=tau, hhat=hhat, r=r, n=n)}

R/hybrid.R

@@ -162,7 +162,7 @@ dAoYda = function(a, vars, pars, FoIpar){with(as.list(c(vars,pars)),{
 solve_dAoYda = function(h, FoIpar, r=1/200, tau=0, Amax=730, dt=1, n=8){
   tms = seq(0, Amax, by = dt)
   prms = c(h=h, r=r, tau=tau, n=n)
-  inits = c(m=0, x=0, y=0)
+  inits = c(m=1e-8, x=0, y=0)
   data.frame(deSolve::ode(inits, times=tms, dAoYda, prms, FoIpar=FoIpar))
 }
 

R/prevalence.R

@@ -3,15 +3,15 @@
 #'
 #' @param a the host cohort age
 #' @param FoIpar parameters that define an FoI function
-#' @param hhat a local scaling parameter for the FoI
 #' @param tau the cohort birthday
+#' @param hhat a local scaling parameter for the FoI
 #' @param r the clearance rate for a simple infection
 #'
 #' @return a [numeric] vector of length(a)
 #' @export
 #'
 truePRa = function(a, FoIpar, tau=0, hhat=1, r=1/200){
-  1-exp(-meanMoI(a, FoIpar, hhat, tau, r))
+  1-exp(-meanMoI(a, FoIpar, tau, hhat, r))
 }
 
 

docs/pkgdown.yml

@@ -9,7 +9,7 @@ articles:
   Parasite_Densities: Parasite_Densities.html
   plotCounts: plotCounts.html
   prevalence: prevalence.html
-last_built: 2024-02-13T22:42Z
+last_built: 2024-02-13T23:35Z
 urls:
   reference: https://dd-harp.github.io/pf.memory/reference
   article: https://dd-harp.github.io/pf.memory/articles

vignettes/Hybrid.Rmd

@@ -194,12 +194,12 @@ Code exists in `pf.memory` to compute prevalence in three ways:
 
 + From the hybrid model `solve_dm` and $p(t)= 1 - exp(-m(t))$
 
-
 ```{r, fig.height=4, fig.width=6, echo=F}
-foiP3 = list(hbar = 5/365, agePar = par_type2Age(), seasonPar = par_sinSeason(), trendPar = par_flatTrend()) 
+foiP3 = list(hbar=1, agePar = par_type2Age(), seasonPar = par_sinSeason(), trendPar = par_flatTrend()) 
 tm = 0:1825
 Xt= solve_dpda(5/365, foiP3, Amax=5*365)
-plot(tm, truePRa(tm, foiP3), type = "l", ylab = expression(p(a)), xlab = "Age (Days)", lwd=5, col = "red")
+tPR = truePRa(tm, foiP3, hhat=5/365)
+plot(tm, tPR, type = "l", ylab = expression(p(a)), xlab = "Age (Days)", lwd=5, col = "red")
 lines(Xt$time, 1-exp(-Xt$m), col = "white", lwd=2)
 lines(Xt$time, Xt$p, col = "blue", lwd=2, lty=2)
 ```
@@ -291,7 +291,7 @@ $$
 To verify, we can compute the moment directly. The function `solve_dAoYda` gives solutions: 
 
 ```{r}
-solve_dAoYda(5/365, foiP3, Amax=5*365, dt=5, n=5) -> mt
+solve_dAoYda(5/365, foiP3, Amax=5*365, dt=5, n=9) -> mt
 ```
 
 The moments can be computed directly using `momentAoY`
@@ -301,9 +301,6 @@ moment1y = momentAoY(aa, foiP3, hhat=5/365)
 ```
 
 
-```{r, echo=F, eval=F}
-devtools::load_all()
-```
 
 The following plots the first moment computed both ways:
 
@@ -312,9 +309,9 @@ plot(aa, moment1y, type = "l", lwd=2, xlab = "Age (in Days)")
 with(mt, lines(time, y, col = "darkred", lwd=2))
 ```
 
-The errors are minor. After 5 years,  the difference between the two methods is less than two days. 
+<!-- The errors are minor. After 5 years,  the difference between the two methods is less than two days.  -->
 
-```{r, fig.height=4, fig.width=6, echo=F}
+```{r, fig.height=4, fig.width=6, echo=F, eval=F}
 plot(moment1y-mt$y[-1], type = "l", xlab = "Age (in Days)", ylab = "Errors")
 ```