bugfix for posterior sampler

Author: Pedro Arruda

Cargo.lock

@@ -1,6 +1,6 @@
 [package]
 name = "ppca_rs"
-version = "0.5.0"
+version = "0.5.1"
 edition = "2021"
 publish = false
 

Cargo.toml

@@ -5,12 +5,12 @@
 
 transform = np.random.binomial(1.0, 0.1, size=(200, 16))
 real_model = PPCAModel(
-    transform=np.matrix(transform, dtype="float64").T.T,
+    transform=np.matrix(transform, dtype="float64"),
     isotropic_noise=0.1,
     mean=np.zeros((200, 1), dtype="float64"),
 )
 
-print("Generating syntetic sample")
+print("Generating synthetic sample")
 sample = real_model.sample(100_000, 0.2)
 
 print("Initializing model")

examples/big_toy_model.py

@@ -1,6 +1,6 @@
 [package]
 name = "ppca"
-version = "0.5.0"
+version = "0.5.1"
 edition = "2021"
 authors = ["Pedro Bittencourt Arruda <[email protected]>"]
 description = "Rust implementation of the Probabilistic Principal Component Analysis model"

ppca/Cargo.toml

@@ -596,11 +596,23 @@ impl Distribution<DVector<f64>> for PosteriorSampler {
     where
         R: Rng + ?Sized,
     {
+        // State noise:
         let standard: DVector<f64> = (0..self.state.len())
             .map(|_| rand_distr::StandardNormal.sample(rng))
             .collect::<Vec<_>>()
             .into();
-        self.model.mean() + self.model.transform() * (&self.state + &self.cholesky_l * standard)
+        // Output noise:
+        let noise: DVector<f64> = (0..self.model.output_size())
+            .map(|_| {
+                let standard: f64 = rand_distr::StandardNormal.sample(rng);
+                self.model.0.output_covariance.isotropic_noise * standard
+            })
+            .collect::<Vec<_>>()
+            .into();
+
+        noise
+            + self.model.mean()
+            + self.model.transform() * (&self.state + &self.cholesky_l * standard)
     }
 }
 

ppca/src/ppca_model.rs

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "ppca-rs"
-version = "0.5.0"
+version = "0.5.1"
 requires-python = ">=3.7"
 description = "Python+Rust implementation of the Probabilistic Principal Component Analysis model"
 readme = "readme.md"

pyproject.toml

@@ -8,7 +8,7 @@
 import numpy as np
 
 
-__version__ = "0.5.0"
+__version__ = "0.5.1"
 
 
 @dataclass(frozen=True)