diff --git a/dev/articles/fitting-dists-with-stan.html b/dev/articles/fitting-dists-with-stan.html index 5a62279..1febd71 100644 --- a/dev/articles/fitting-dists-with-stan.html +++ b/dev/articles/fitting-dists-with-stan.html @@ -151,31 +151,35 @@

pwindow = pwindows, swindow = swindows, obs_time = obs_times, - observed_delay = samples + observed_delay = samples, + observed_delay_upper = samples + swindows ) head(delay_data) -
##   pwindow swindow obs_time observed_delay
-## 1       1       1       10              2
-## 2       1       1       10              1
-## 3       1       1       10              8
-## 4       1       1       10              6
-## 5       1       1       10              5
-## 6       1       1       10              4
+
##   pwindow swindow obs_time observed_delay observed_delay_upper
+## 1       1       1       10              2                    3
+## 2       1       1       10              1                    2
+## 3       1       1       10              8                    9
+## 4       1       1       10              6                    7
+## 5       1       1       10              5                    6
+## 6       1       1       10              4                    5
 # Aggregate to unique combinations and count occurrences
 # Aggregate to unique combinations and count occurrences
 delay_counts <- delay_data |>
-  summarise(n = n(), .by = c(pwindow, swindow, obs_time, observed_delay))
+  summarise(
+    n = n(),
+    .by = c(pwindow, swindow, obs_time, observed_delay, observed_delay_upper)
+  )
 
 head(delay_counts)
-
##   pwindow swindow obs_time observed_delay   n
-## 1       1       1       10              2 160
-## 2       1       1       10              1  96
-## 3       1       1       10              8  59
-## 4       1       1       10              6 113
-## 5       1       1       10              5 119
-## 6       1       1       10              4 150
+
##   pwindow swindow obs_time observed_delay observed_delay_upper   n
+## 1       1       1       10              2                    3 160
+## 2       1       1       10              1                    2  96
+## 3       1       1       10              8                    9  59
+## 4       1       1       10              6                    7 113
+## 5       1       1       10              5                    6 119
+## 6       1       1       10              4                    5 150
 # Compare the samples with and without secondary censoring to the true
 # distribution
@@ -265,8 +269,8 @@ 

), chains = 4, parallel_chains = 4, - show_messages = FALSE, - refresh = 0 + refresh = ifelse(interactive(), 50, 0), + show_messages = interactive() ) naive_fit

## Warning: NAs introduced by coercion
@@ -295,9 +299,9 @@

data { int<lower=0> N; // number of observations array[N] int<lower=0> y; // observed delays + array[N] int<lower=0> y_upper; // observed delays upper bound array[N] int<lower=0> n; // number of occurrences for each delay array[N] int<lower=0> pwindow; // primary censoring window - array[N] int<lower=0> swindow; // secondary censoring window array[N] int<lower=0> D; // maximum delay } transformed data { @@ -316,7 +320,7 @@

for (i in 1:N) { target += n[i] * primary_censored_dist_lpmf( y[i] | 1, {mu, sigma}, - pwindow[i], swindow[i], D[i], + pwindow[i], y_upper[i], D[i], 1, primary_params ); } @@ -334,28 +338,22 @@

primarycensoreddist_fit <- primarycensoreddist_model$sample( data = list( y = delay_counts$observed_delay, + y_upper = delay_counts$observed_delay_upper, n = delay_counts$n, pwindow = delay_counts$pwindow, - swindow = delay_counts$swindow, D = delay_counts$obs_time, N = nrow(delay_counts) ), chains = 4, parallel_chains = 4, - init = list( # we use this to resolve initialisation issues - list(mu = 1.5, sigma = 0.6), - list(mu = 1.5, sigma = 0.4), - list(mu = 1.5, sigma = 0.3), - list(mu = 1.5, sigma = 0.55) - ), - refresh = 0, - show_messages = FALSE + refresh = ifelse(interactive(), 50, 0), + show_messages = interactive() ) primarycensoreddist_fit
##  variable     mean   median   sd  mad       q5      q95 rhat ess_bulk ess_tail
-##     lp__  -2134.42 -2134.13 0.97 0.70 -2136.43 -2133.49 1.00      732      952
-##     mu        1.53     1.53 0.05 0.05     1.46     1.61 1.01      620      686
-##     sigma     0.77     0.77 0.04 0.04     0.71     0.83 1.00      633      641
+## lp__ -2134.45 -2134.16 0.97 0.71 -2136.44 -2133.50 1.00 1741 2206 +## mu 1.53 1.53 0.05 0.05 1.46 1.61 1.00 1386 1540 +## sigma 0.77 0.77 0.04 0.04 0.71 0.83 1.00 1373 1786

We see that the model has converged and the diagnostics look good. We also see that the posterior means are very near the true parameters and the 90% credible intervals include the true parameters.