From a1e6eb85a093b14139340581d528ae8f381cc718 Mon Sep 17 00:00:00 2001 From: Emil Koutanov Date: Thu, 4 Jan 2024 17:12:15 +1100 Subject: [PATCH] Added Beta and Pareto --- README.md | 2 +- examples/sample_beta.rs | 35 +++++++++++++++++++++++++++++++++++ examples/sample_pareto.rs | 35 +++++++++++++++++++++++++++++++++++ justfile | 16 ++++++++++++---- src/beta.rs | 21 +++++++++++++++++++++ src/lib.rs | 2 ++ src/pareto.rs | 21 +++++++++++++++++++++ 7 files changed, 127 insertions(+), 5 deletions(-) create mode 100644 examples/sample_beta.rs create mode 100644 examples/sample_pareto.rs create mode 100644 src/beta.rs create mode 100644 src/pareto.rs diff --git a/README.md b/README.md index caebcd6..8b6d3c6 100644 --- a/README.md +++ b/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. diff --git a/examples/sample_beta.rs b/examples/sample_beta.rs new file mode 100644 index 0000000..6c9b572 --- /dev/null +++ b/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 = 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); + } +} \ No newline at end of file diff --git a/examples/sample_pareto.rs b/examples/sample_pareto.rs new file mode 100644 index 0000000..0c25066 --- /dev/null +++ b/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 = 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); + } +} \ No newline at end of file diff --git a/justfile b/justfile index 9b9aa2e..662d098 100644 --- a/justfile +++ b/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' diff --git a/src/beta.rs b/src/beta.rs new file mode 100644 index 0000000..d5f59a0 --- /dev/null +++ b/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) + } +} \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index c92ec1b..8e5941d 100644 --- a/src/lib.rs +++ b/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 { diff --git a/src/pareto.rs b/src/pareto.rs new file mode 100644 index 0000000..c93c896 --- /dev/null +++ b/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) + } +} \ No newline at end of file