Added Beta and Pareto

README.md

@@ -6,7 +6,7 @@ Markov chain Monte Carlo sampling using the _Independence Metropolis-Hastings_ a
 [![docs.rs](https://img.shields.io/badge/docs.rs-metromc-blue?style=flat-square&logo=docs.rs)](https://docs.rs/metromc)
 [![Build Status](https://img.shields.io/github/actions/workflow/status/obsidiandynamics/metromc/master.yml?branch=master&style=flat-square&logo=github)](https://github.com/obsidiandynamics/metromc/actions/workflows/master.yml)
 
-Uses the [tinyrand](https://github.com/obsidiandynamics/tinyrand) RNG to sample at rate of ~50M samples/sec.
+Uses the [tinyrand](https://github.com/obsidiandynamics/tinyrand) RNG to sample at a rate of ~50M samples/sec.
 
 # Example
 Draw samples from the Gaussian distribution.

examples/sample_beta.rs

@@ -0,0 +1,35 @@
+use std::ops::RangeInclusive;
+
+use tinyrand::Wyrand;
+use metromc::beta::Beta;
+
+use metromc::sampler::{Config, Sampler};
+
+fn main() {
+    const ALPHA: f64 = 0.5;
+    const BETA: f64 = 0.5;
+    const RANGE: RangeInclusive<f64> = 0.0..=1.0;
+    const NUM_BUCKETS: usize = 1_000;
+
+    let sampler = Sampler::new(Config {
+        rand: Wyrand::default(),
+        dist: Beta::new(ALPHA, BETA),
+        range: RANGE,
+    });
+
+    let mut buckets = [0_usize; NUM_BUCKETS];
+    let span = RANGE.end() - RANGE.start();
+    for sample in sampler.skip(100).take(10_000_000) {
+        // println!("{rand:.6}");
+        let bucket = ((sample - RANGE.start()) / span * NUM_BUCKETS as f64) as usize;
+        buckets[bucket] += 1;
+    }
+
+    println!("bucket start, count");
+    println!("------------, -----");
+
+    for (bucket, count) in buckets.iter().enumerate() {
+        let bucket_start = RANGE.start() + bucket as f64 / NUM_BUCKETS as f64 * span;
+        println!("{bucket_start:.3}, {}", count);
+    }
+}

examples/sample_pareto.rs

@@ -0,0 +1,35 @@
+use std::ops::RangeInclusive;
+
+use tinyrand::Wyrand;
+
+use metromc::pareto::Pareto;
+use metromc::sampler::{Config, Sampler};
+
+fn main() {
+    const SHAPE: f64 = 2.0;
+    const RATE: f64 = 1.0;
+    const RANGE: RangeInclusive<f64> = 0.0..=16.0;
+    const NUM_BUCKETS: usize = 1_000;
+
+    let sampler = Sampler::new(Config {
+        rand: Wyrand::default(),
+        dist: Pareto::new(SHAPE, RATE),
+        range: RANGE,
+    });
+
+    let mut buckets = [0_usize; NUM_BUCKETS];
+    let span = RANGE.end() - RANGE.start();
+    for sample in sampler.skip(100).take(10_000_000) {
+        // println!("{rand:.6}");
+        let bucket = ((sample - RANGE.start()) / span * NUM_BUCKETS as f64) as usize;
+        buckets[bucket] += 1;
+    }
+
+    println!("bucket start, count");
+    println!("------------, -----");
+
+    for (bucket, count) in buckets.iter().enumerate() {
+        let bucket_start = RANGE.start() + bucket as f64 / NUM_BUCKETS as f64 * span;
+        println!("{bucket_start:.3}, {}", count);
+    }
+}

justfile

@@ -1,22 +1,30 @@
 _help:
     @just --list
 
-# list values of a Gaussian pdf
+# list values of the Gaussian pdf
 gaussian_pdf:
     cargo run --release --example gaussian_pdf
 
-# list values of a Gamma pdf
+# list values of the Gamma pdf
 gamma_pdf:
     cargo run --release --example gamma_pdf
 
-# sample from a Gaussian distribution
+# sample from the Gaussian distribution
 sample_gaussian:
     cargo run --release --example sample_gaussian
 
-# sample from a Gamma distribution
+# sample from the Gamma distribution
 sample_gamma:
     cargo run --release --example sample_gamma
 
+# sample from the Pareto distribution
+sample_pareto:
+    cargo run --release --example sample_pareto
+
+# sample from the Beta distribution
+sample_beta:
+    cargo run --release --example sample_beta
+
 # run Criterion bechmarks
 bench:
     bash -c 'type cargo-criterion >/dev/null 2>&1 || cargo install cargo-criterion'

src/beta.rs

@@ -0,0 +1,21 @@
+use statrs::distribution::Continuous;
+use crate::Pdf;
+
+pub struct Beta {
+    dist: statrs::distribution::Beta,
+}
+impl Beta {
+    #[inline(always)]
+    pub fn new(alpha: f64, beta: f64) -> Self {
+        Self {
+            dist: statrs::distribution::Beta::new(alpha, beta).unwrap(),
+        }
+    }
+}
+
+impl Pdf for Beta {
+    #[inline(always)]
+    fn prob(&self, x: f64) -> f64 {
+        self.dist.pdf(x)
+    }
+}

src/lib.rs

@@ -1,5 +1,7 @@
+pub mod beta;
 pub mod gamma;
 pub mod gaussian;
+pub mod pareto;
 pub mod sampler;
 
 pub trait Pdf {

src/pareto.rs

@@ -0,0 +1,21 @@
+use statrs::distribution::Continuous;
+use crate::Pdf;
+
+pub struct Pareto {
+    dist: statrs::distribution::Pareto,
+}
+impl Pareto {
+    #[inline(always)]
+    pub fn new(shape: f64, rate: f64) -> Self {
+        Self {
+            dist: statrs::distribution::Pareto::new(shape, rate).unwrap(),
+        }
+    }
+}
+
+impl Pdf for Pareto {
+    #[inline(always)]
+    fn prob(&self, x: f64) -> f64 {
+        self.dist.pdf(x)
+    }
+}