1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include lib.stan
# This is a deprecated way to write a comment!
functions {
/** Solution to drug concentration ODE given inital concentration, time elapsed,
* and diffusion parameters.
*/
vector[] drug_conc(vector y_cen, vector y_per, vector kappa_cen,
vector kappa_per, vector ts, int N) {
vector[N] y_cen_out = exp(-kappa_cen .* ts) .* y_cen;
vector[N] y_per_out = (kappa_cen ./ (kappa_per - kappa_cen))
.* (exp(-kappa_cen .* ts) - exp(-kappa_per .* ts))
.* y_cen + exp(-kappa_per .* ts) .* y_per;
return { y_cen_out, y_per_out };
}
// Functor with appropriate signature for algebraic solver.
vector f(vector y, vector kappas, real[] x_r, int[] x_i) {
// Unpack x_r.
real delta = x_r[1];
real tau = x_r[2];
// Unpack x_i.
int n = x_i[1];
// All of the intervals are tau.
vector[n] ts = rep_vector(tau, n);
/* The first n entries of y are the concentrations in the central chambers,
* while the last n entries of y are the concentrations in the peripheral
* chambers. Likewise for kappa.
*/
vector[n] y_cen = y[:n];
vector[n] y_per = y[(n+1):];
vector[n] kappa_cen = kappas[:n];
vector[n] kappa_per = kappas[(n+1):];
// Calculate the concentrations after tau.
vector[n] y_res[2] = drug_conc(y_cen, y_per, kappa_cen, kappa_per, ts, n);
/* If a steady state has been reached, the difference between the
* concentration after tau (with a dose delta) should be the same as the
* current state.
*/
return append_row(y_res[1] + rep_vector(delta, n), y_res[2]) - y;
}
#include addl_functions.stanfunctions // With a comment!
}
data {
real<lower=0> delta; // Dosage
real<lower=0> tau; // Dose interval
int n; // Number of patients
int m; // Number of observations
vector<lower=0,upper=tau>[m] ts; // Times of observations (since last dose)
int<lower=0,upper=n> idx[m]; // Patient corresponding to observation
vector<lower=0>[m] obs; // Observed concetrations
vector<lower=0>[n] y_guess_cen; // Guess for central chamber conc.
vector<lower=0>[n] y_guess_per; // Guess for peripheral chamber conc.
}
transformed data {
// Reshape the data for the algebra solver.
vector[2*n] y_guess = append_row(y_guess_cen, y_guess_per);
real x_r[2] = { delta, tau };
int x_i[1] = { n };
}
parameters {
vector<lower=0>[n] kappa_cen; // Dispersion from central chamber
vector<lower=0>[n] kappa_per; // Dispersion from peripheral chamber
array[n] real dummy_variables; // Dummy variables
complex complex_dummy_variable; // complex dummy variable
}
transformed parameters {
// Reshape dispersion parameters for algebra solver.
vector[2*n] kappas = append_row(kappa_cen, kappa_per);
// Get the steady-state for each patient.
vector[2*n] y_steady = algebra_solver(f, y_guess, kappas, x_r, x_i);
vector[n] y_steady_cen = y_steady[:n];
vector[n] y_steady_per = y_steady[(n+1):];
/* Get the concentrations in the central chamber at each time observation,
* given the currently sampled diffusion parameters.
*/
vector[m] y_true[2] = drug_conc(y_steady_cen[idx], y_steady_per[idx],
kappa_cen[idx], kappa_per[idx], ts, m);
// Add a constant to the complex dummy variable
complex_dummy_variable += 3i - 40e-3i + 1e10i
}
model {
target += normal_lpdf(kappa_cen | get_real(0 + 0.0i), 1.0/4) T[0, ];
target += normal_lpdf(kappa_per | get_imag(0 + 0i), pi()/4) T[0, ];
obs ~ lognormal(log(y_true[2]), 1.0/4);
}
generated quantities {
real fake_obs[m] = lognormal_rng(log(y_true[2]), 1.0/4);
}
|
