pharmsol 0.27.1

Rust library for solving analytic and ode-defined pharmacometric models.
Documentation 
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

use pharmsol::data::error_model::AssayErrorModel;
use pharmsol::*;

/// Test the particle filter (SDE) likelihood estimation
/// This test runs the particle filter multiple times to assess consistency
#[test]
fn test_particle_filter_likelihood() {
    let subject = data::Subject::builder("id1")
        .bolus(0.0, 20.0, "dose")
        .observation(0.2, 16.6434, "cp")
        .observation(0.4, 14.3233, "cp")
        .observation(0.6, 9.8468, "cp")
        .observation(0.8, 9.4177, "cp")
        .observation(1.0, 7.5170, "cp")
        .build();

    let sde = equation::SDE::new(
        |x, p, _t, dx, _rateiv, _cov| {
            dx[0] = -x[0] * x[1]; // ke *x[0]
            dx[1] = -x[1] + p[0]; // mean reverting
        },
        |_p, d| {
            d[0] = 1.0;
            d[1] = 0.01;
        },
        |_p, _t, _cov| lag! {},
        |_p, _t, _cov| fa! {},
        |_p, _t, _cov, x| x[1] = 1.0,
        |x, _p, _t, _cov, y| {
            y[0] = x[0];
        },
        10000,
    )
    .with_nstates(2)
    .with_ndrugs(1)
    .with_nout(1);

    let sde = sde
        .with_metadata(
            equation::metadata::new("particle_filter_test")
                .kind(equation::ModelKind::Sde)
                .parameters(["ke0"])
                .states(["central", "ke_latent"])
                .outputs(["cp"])
                .route(
                    equation::Route::bolus("dose")
                        .to_state("central")
                        .inject_input_to_destination(),
                )
                .particles(10000),
        )
        .expect("particle filter metadata should validate");

    let ems = AssayErrorModels::default()
        .add(
            0,
            AssayErrorModel::additive(ErrorPoly::new(0.5, 0.0, 0.0, 0.0), 0.0),
        )
        .unwrap();

    // Run the particle filter multiple times to get statistics
    const NUM_RUNS: usize = 10;
    let mut likelihoods = Vec::with_capacity(NUM_RUNS);

    for i in 0..NUM_RUNS {
        let parameters = Parameters::with_model(&sde, [("ke0", 1.0)])
            .expect("particle filter parameters should validate");
        let ll = sde
            .estimate_log_likelihood(&subject, &parameters, &ems)
            .unwrap()
            .exp();
        println!("Run {}: likelihood = {}", i + 1, ll);
        likelihoods.push(ll);
    }

    // Calculate mean and standard deviation
    let mean: f64 = likelihoods.iter().sum::<f64>() / NUM_RUNS as f64;
    let variance: f64 =
        likelihoods.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / NUM_RUNS as f64;
    let std_dev = variance.sqrt();

    println!("\n=== Particle Filter Likelihood Statistics ===");
    println!("Number of runs: {}", NUM_RUNS);
    println!("Likelihoods: {:?}", likelihoods);
    println!("Mean likelihood: {}", mean);
    println!("Std deviation: {}", std_dev);
    println!(
        "Min: {}",
        likelihoods.iter().cloned().fold(f64::INFINITY, f64::min)
    );
    println!(
        "Max: {}",
        likelihoods
            .iter()
            .cloned()
            .fold(f64::NEG_INFINITY, f64::max)
    );

    // The likelihood should be a reasonable negative value (log-likelihood)
    // This assertion is loose since particle filters have stochastic variation
    assert!(mean.is_finite(), "Mean likelihood should be finite");
}