rhmm 0.0.2

A rust implementation of hidden markov 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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

//! Sampling utilities

use crate::errors::{HmmError, Result};
use ndarray::Array1;
use rand::Rng;

/// Sample from a discrete distribution
///
/// # Arguments
///
/// * `probs` - Probability distribution
/// * `rng` - Random number generator
///
/// # Returns
///
/// Sampled index
pub fn sample_discrete<R: Rng>(probs: &Array1<f64>, rng: &mut R) -> Result<usize> {
    let sum: f64 = probs.sum();
    if (sum - 1.0).abs() > 1e-6 {
        return Err(HmmError::InvalidProbability(format!(
            "Probabilities must sum to 1.0, got {}",
            sum
        )));
    }

    let mut cumsum = 0.0;
    let rand_val: f64 = rng.random();

    for (i, &p) in probs.iter().enumerate() {
        cumsum += p;
        if rand_val <= cumsum {
            return Ok(i);
        }
    }

    // Fallback to last index (handles floating point errors)
    Ok(probs.len() - 1)
}

/// Sample from a multivariate normal distribution
///
/// # Arguments
///
/// * `mean` - Mean vector
/// * `covar` - Covariance (diagonal elements only for now)
/// * `rng` - Random number generator
///
/// # Returns
///
/// Sampled vector
pub fn sample_gaussian<R: Rng>(
    mean: &Array1<f64>,
    covar: &Array1<f64>,
    rng: &mut R,
) -> Result<Array1<f64>> {
    use rand_distr::{Distribution, Normal};

    let n_features = mean.len();
    let mut sample = Array1::zeros(n_features);

    for i in 0..n_features {
        let std_dev = covar[i].sqrt();
        let normal = Normal::new(mean[i], std_dev).map_err(|e| {
            HmmError::InvalidParameter(format!("Invalid normal distribution parameters: {}", e))
        })?;
        sample[i] = normal.sample(rng);
    }

    Ok(sample)
}

#[cfg(test)]
mod tests {
    use super::*;
    use ndarray::array;
    use rand::rngs::StdRng;
    use rand::SeedableRng;

    #[test]
    fn test_sample_discrete() {
        let mut rng = StdRng::seed_from_u64(42);
        let probs = array![0.5, 0.3, 0.2];

        let mut counts = vec![0; 3];
        for _ in 0..1000 {
            let idx = sample_discrete(&probs, &mut rng).unwrap();
            counts[idx] += 1;
        }

        // Check that sampling roughly follows the distribution
        assert!(counts[0] > counts[1]);
        assert!(counts[1] > counts[2]);
    }

    #[test]
    fn test_sample_discrete_invalid_sum() {
        let mut rng = StdRng::seed_from_u64(42);
        let probs = array![0.5, 0.3, 0.3];

        assert!(sample_discrete(&probs, &mut rng).is_err());
    }

    #[test]
    fn test_sample_gaussian() {
        let mut rng = StdRng::seed_from_u64(42);
        let mean = array![0.0, 1.0];
        let covar = array![1.0, 0.5];

        let sample = sample_gaussian(&mean, &covar, &mut rng).unwrap();
        assert_eq!(sample.len(), 2);
    }

    #[test]
    fn test_sample_gaussian_multiple() {
        let mut rng = StdRng::seed_from_u64(42);
        let mean = array![5.0];
        let covar = array![1.0];

        let mut samples = Vec::new();
        for _ in 0..100 {
            let sample = sample_gaussian(&mean, &covar, &mut rng).unwrap();
            samples.push(sample[0]);
        }

        let sample_mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
        // Mean should be roughly around 5.0
        assert!((sample_mean - 5.0).abs() < 1.0);
    }
}