feat: add binomial distribution calc #1

Merged
pitbuster merged 5 commits from binom into main 2025-12-08 15:54:10 -03:00
4 changed files with 316 additions and 115 deletions
Showing only changes of commit 86f822bf34 - Show all commits

View file

@ -10,4 +10,4 @@ hyper-description = Hypergeometric distribution measures the probability of gett
success-probability = Success probability success-probability = Success probability
trials-number = Number of trials trials-number = Number of trials
successes-number = Number of successes successes-number = Number of successes
binom-description = . binom-description = Binomial distribution measures the probability of getting a given amount of successes in a sequence of experiments. For example, if you flip 5 (number of trials = 5) balanced coins (success probability = 0.5), the distribution describes the probability of having a given number of heads (successes number = X).

View file

@ -10,4 +10,4 @@ hyper-description = La distribución hipergeométrica mide la probabilidad de ob
success-probability = Probabilidad de éxito success-probability = Probabilidad de éxito
trials-number = Número de intentos trials-number = Número de intentos
successes-number = Número de éxitos successes-number = Número de éxitos
binom-description = . binom-description = La distribución binomial mide la probabilidad de obtener un cierto número de éxitos en una secuencia de experimentos. Por ejemplo, si lanzas 5 (número de intentos = 5) monedas justas (probabilidad de éxito = 0.5), la distribución describe la probabilidad de obtener un cierto número de monedas (número de éxitos = X).

View file

@ -2,79 +2,179 @@
use std::{collections::HashMap, iter::repeat}; use std::{collections::HashMap, iter::repeat};
#[derive(Default)] #[derive(Debug)]
pub struct HyperGeometricInput {
population_size: u8,
successes: u8,
sample_size: u8,
sample_successes: u8,
}
impl HyperGeometricInput {
pub fn new(
population_size: u8,
successes: u8,
sample_size: u8,
sample_successes: u8,
) -> Option<Self> {
if successes > population_size
|| sample_size > population_size
|| sample_successes > sample_size
{
None
} else {
Some(Self {
population_size,
successes,
sample_size,
sample_successes,
})
}
}
}
/// Result of hypergeometric probability calculation.
#[derive(Default, Debug, PartialEq)]
pub struct HyperGeometricProb { pub struct HyperGeometricProb {
/// Probability of getting exactly X successes in the sample.
pub exactly: f64, pub exactly: f64,
/// Probability of getting strictly less than X successes in the sample.
pub less_than: f64, pub less_than: f64,
/// Probability of getting less than or exactly X successes in the sample.
pub less_or_equal: f64, pub less_or_equal: f64,
/// Probability of getting strictly more X successes in the sample.
pub greater_than: f64, pub greater_than: f64,
/// Probability of getting more than or exactly X successes in the sample.
pub greater_or_equal: f64, pub greater_or_equal: f64,
} }
pub fn hyper_geometric( pub fn hyper_geometric(input: HyperGeometricInput) -> HyperGeometricProb {
population_size: u8, let exactly = hyper_geometric_exactly(&input);
successes: u8, let (less_than, less_or_equal, greater_or_equal, greater_than) = if input.sample_successes
sample_size: u8, < input.sample_size / 2
sample_successes: u8,
) -> Option<HyperGeometricProb> {
if successes > population_size
|| sample_size > population_size
|| sample_successes > sample_size
{ {
None let less_than = (0..input.sample_successes)
.map(|i| {
hyper_geometric_exactly(&HyperGeometricInput {
population_size: input.population_size,
successes: input.successes,
sample_size: input.sample_size,
sample_successes: i,
})
})
.sum::<f64>()
.abs();
let less_or_equal = less_than + exactly;
let greater_or_equal = (1.0 - less_than).abs();
let greater_than = (1.0 - less_or_equal).abs();
(less_than, less_or_equal, greater_or_equal, greater_than)
} else { } else {
let exactly = let greater_than = (input.sample_successes + 1..=input.sample_size.min(input.successes))
hyper_geometric_exactly(population_size, successes, sample_size, sample_successes); .map(|i| {
let (less_than, less_or_equal, greater_or_equal, greater_than) = hyper_geometric_exactly(&HyperGeometricInput {
if sample_successes < sample_size / 2 { population_size: input.population_size,
let less_than = (0..sample_successes) successes: input.successes,
.map(|i| hyper_geometric_exactly(population_size, successes, sample_size, i)) sample_size: input.sample_size,
.sum::<f64>() sample_successes: i,
.abs(); })
let less_or_equal = less_than + exactly; })
let greater_or_equal = (1.0 - less_than).abs(); .sum::<f64>()
let greater_than = (1.0 - less_or_equal).abs(); .abs();
(less_than, less_or_equal, greater_or_equal, greater_than) let greater_or_equal = greater_than + exactly;
} else { let less_or_equal = (1.0 - greater_than).abs();
let greater_than = (sample_successes + 1..=sample_size) let less_than = (1.0 - greater_or_equal).abs();
.map(|i| hyper_geometric_exactly(population_size, successes, sample_size, i)) (less_than, less_or_equal, greater_or_equal, greater_than)
.sum::<f64>() };
.abs(); HyperGeometricProb {
let greater_or_equal = greater_than + exactly; exactly,
let less_or_equal = (1.0 - greater_than).abs(); less_than,
let less_than = (1.0 - greater_or_equal).abs(); less_or_equal,
(less_than, less_or_equal, greater_or_equal, greater_than) greater_than,
}; greater_or_equal,
Some(HyperGeometricProb {
exactly,
less_than,
less_or_equal,
greater_than,
greater_or_equal,
})
} }
} }
#[derive(Default)] pub struct BinomialInput {
pub struct BinomialProb {
pub exactly: f64,
pub less_than: f64,
pub less_or_equal: f64,
pub greater_than: f64,
pub greater_or_equal: f64,
}
pub fn binomial(
success_probability: f64, success_probability: f64,
trials_number: u8, trials_number: u8,
successes_number: u8, successes_number: u8,
) -> Option<BinomialProb> { }
if successes_number > trials_number {
None impl BinomialInput {
} else { pub fn new(success_probability: f64, trials_number: u8, successes_number: u8) -> Option<Self> {
let (p_powers, pc_powers) = powers(success_probability, trials_number); if success_probability < 0.0
let exactly = binom_exactly(success_probability, trials_number, successes_number); || success_probability > 1.0
None || successes_number > trials_number
{
None
} else {
Some(Self {
success_probability,
trials_number,
successes_number,
})
}
}
}
#[derive(Default, Debug, PartialEq)]
pub struct BinomialProb {
pub exactly: f64,
pub less_than: f64,
pub less_or_equal: f64,
pub greater_than: f64,
pub greater_or_equal: f64,
}
pub fn binomial(input: BinomialInput) -> BinomialProb {
let (p_powers, pc_powers) = powers(input.success_probability, input.trials_number);
let exactly = binomial_exactly(&input, &p_powers, &pc_powers);
let (less_than, less_or_equal, greater_or_equal, greater_than) =
if input.successes_number < input.trials_number / 2 {
let less_than = (0..input.successes_number)
.map(|i| {
binomial_exactly(
&BinomialInput {
success_probability: input.success_probability,
trials_number: input.trials_number,
successes_number: i,
},
&p_powers,
&pc_powers,
)
})
.sum::<f64>()
.abs();
let less_or_equal = less_than + exactly;
let greater_or_equal = (1.0 - less_than).abs();
let greater_than = (1.0 - less_or_equal).abs();
(less_than, less_or_equal, greater_or_equal, greater_than)
} else {
let greater_than = (input.successes_number + 1..=input.trials_number)
.map(|i| {
binomial_exactly(
&BinomialInput {
success_probability: input.success_probability,
trials_number: input.trials_number,
successes_number: i,
},
&p_powers,
&pc_powers,
)
})
.sum::<f64>()
.abs();
let greater_or_equal = greater_than + exactly;
let less_or_equal = (1.0 - greater_than).abs();
let less_than = (1.0 - greater_or_equal).abs();
(less_than, less_or_equal, greater_or_equal, greater_than)
};
BinomialProb {
exactly,
less_than,
less_or_equal,
greater_than,
greater_or_equal,
} }
} }
@ -83,65 +183,49 @@ pub fn binomial(
/// ///
/// The formula is choose(successes, sample_successes) * choose(population_size - successes, /// The formula is choose(successes, sample_successes) * choose(population_size - successes,
/// sample_size - sample_successes) / choose(population_size, sample_size) /// sample_size - sample_successes) / choose(population_size, sample_size)
fn hyper_geometric_exactly( fn hyper_geometric_exactly(input: &HyperGeometricInput) -> f64 {
population_size: u8, if input.population_size == input.successes {
successes: u8, return if input.sample_successes == input.sample_size {
sample_size: u8,
sample_successes: u8,
) -> f64 {
if population_size == successes {
return if sample_successes == sample_size {
1.0 1.0
} else { } else {
0.0 0.0
}; };
} }
if successes == 0 { if input.successes == 0 {
return if sample_successes == 0 { 1.0 } else { 0.0 }; return if input.sample_successes == 0 {
1.0
} else {
0.0
};
} }
// On top we have: successes!, (population_size - successes)!, sample_size! and // On top we have: successes!, (population_size - successes)!, sample_size! and
// (population_size - sample_size)! // (population_size - sample_size)!
let top_factors = (1..=successes) let top_factors = (1..=input.successes)
.chain(1..=(population_size - successes)) .chain(1..=(input.population_size - input.successes))
.chain(1..=sample_size) .chain(1..=input.sample_size)
.chain(1..=(population_size - sample_size)) .chain(1..=(input.population_size - input.sample_size))
.flat_map(|n| factorize(n)) .flat_map(|n| factorize(n))
.fold(HashMap::<u8, u8>::new(), |mut counts, i| { .fold(HashMap::<u8, u8>::new(), group_factors);
*counts.entry(i).or_default() += 1;
counts
});
// On bottom we have: sample_successes!, (successes - sample_successes)! // On bottom we have: sample_successes!, (successes - sample_successes)!
// (sample_size - sample_successes)!, (population_size - successes - sample_size + sample_successes)! // (sample_size - sample_successes)!, (population_size - successes - sample_size + sample_successes)!
// and population_size! // and population_size!
let bot_factors = (1..=sample_successes) let bot_factors = (1..=input.sample_successes)
.chain(1..=(successes - sample_successes)) .chain(1..=(input.successes - input.sample_successes))
.chain(1..=(sample_size - sample_successes)) .chain(1..=(input.sample_size - input.sample_successes))
.chain( .chain(
1..=((population_size as u16 + sample_successes as u16 1..=((input.population_size as u16 + input.sample_successes as u16
- successes as u16 - input.successes as u16
- sample_size as u16) as u8), - input.sample_size as u16) as u8),
) )
.chain(1..=population_size) .chain(1..=input.population_size)
.flat_map(|n| factorize(n)) .flat_map(|n| factorize(n))
.fold(HashMap::<u8, u8>::new(), |mut counts, i| { .fold(HashMap::<u8, u8>::new(), group_factors);
counts.entry(i).and_modify(|count| *count += 1).or_insert(1);
counts
});
let (top_factors, bot_factors) = simplify(top_factors, bot_factors); let (top_factors, bot_factors) = simplify(top_factors, bot_factors);
let top_product: f64 = top_factors let top_product = product(top_factors);
.into_iter() let bot_product = product(bot_factors);
.flat_map(|(f, count)| repeat(f).take(count as usize))
.map(|f| f as f64)
.product();
let bot_product: f64 = bot_factors
.into_iter()
.flat_map(|(f, count)| repeat(f).take(count as usize))
.map(|f| f as f64)
.product();
top_product / bot_product top_product / bot_product
} }
@ -149,19 +233,55 @@ fn hyper_geometric_exactly(
/// Computes the probability of getting exactly `successes_number` within `trials_number` given /// Computes the probability of getting exactly `successes_number` within `trials_number` given
/// that the success probability is `success_probability`. /// that the success probability is `success_probability`.
/// ///
/// The formula is choose(successes_number, trials_number) * (success_probability)^successes_number /// The formula is choose(trials_number, successes_number) * (success_probability)^successes_number
/// * (1 - success_probability)^(trials_number - successes_number) /// * (1 - success_probability)^(trials_number - successes_number)
fn binom_exactly(success_probability: f64, trials_number: u8, successes_number: u8) -> f64 { fn binomial_exactly(input: &BinomialInput, p_powers: &[f64], pc_powers: &[f64]) -> f64 {
0.0 if input.success_probability == 0.0 {
return if input.successes_number == 0 {
1.0
} else {
0.0
};
}
if input.success_probability == 1.0 {
return if input.successes_number == input.trials_number {
1.0
} else {
0.0
};
}
choose(input.trials_number, input.successes_number)
* p_powers[input.successes_number as usize]
* pc_powers[(input.trials_number - input.successes_number) as usize]
} }
fn powers(p: f64, N: u8) -> (Vec<f64>, Vec<f64>) { fn choose(n: u8, k: u8) -> f64 {
let mut p_powers = Vec::with_capacity((N + 1) as usize); // On top we have: n!
let mut pc_powers = Vec::with_capacity((N + 1) as usize); let top_factors = (1..=n)
.flat_map(|n| factorize(n))
.fold(HashMap::<u8, u8>::new(), group_factors);
// On bottom we have: k!, (n - k)!
let bot_factors = (1..=k)
.chain(1..=(n - k))
.flat_map(|n| factorize(n))
.fold(HashMap::<u8, u8>::new(), group_factors);
let (top_factors, bot_factors) = simplify(top_factors, bot_factors);
let top_product = product(top_factors);
let bot_product = product(bot_factors);
top_product / bot_product
}
fn powers(p: f64, n: u8) -> (Vec<f64>, Vec<f64>) {
let mut p_powers = Vec::with_capacity((n + 1) as usize);
let mut pc_powers = Vec::with_capacity((n + 1) as usize);
let mut p_power = 1.0; let mut p_power = 1.0;
let mut pc_power = 1.0; let mut pc_power = 1.0;
for _ in 0..N + 1 { for _ in 0..n + 1 {
p_powers.push(p_power); p_powers.push(p_power);
pc_powers.push(pc_power); pc_powers.push(pc_power);
p_power = p_power * p; p_power = p_power * p;
@ -170,6 +290,19 @@ fn powers(p: f64, N: u8) -> (Vec<f64>, Vec<f64>) {
(p_powers, pc_powers) (p_powers, pc_powers)
} }
fn group_factors(mut counts: HashMap<u8, u8>, i: u8) -> HashMap<u8, u8> {
*counts.entry(i).or_default() += 1;
counts
}
fn product(factors: HashMap<u8, u8>) -> f64 {
factors
.into_iter()
.flat_map(|(f, count)| repeat(f).take(count as usize))
.map(|f| f as f64)
.product()
}
/// Simplify factors for a fraction. /// Simplify factors for a fraction.
/// ///
/// This assumes factors are already prime factors. /// This assumes factors are already prime factors.
@ -259,7 +392,10 @@ fn factorize(n: u8) -> FactorIter<'static> {
#[cfg(test)] #[cfg(test)]
mod test { mod test {
use crate::calc::hyper_geometric_exactly; use crate::calc::{
BinomialInput, BinomialProb, HyperGeometricInput, HyperGeometricProb, binomial,
binomial_exactly, hyper_geometric, hyper_geometric_exactly, powers,
};
use super::factorize; use super::factorize;
@ -280,18 +416,80 @@ mod test {
#[test] #[test]
fn test_hypergeometric_exact_all_successes() { fn test_hypergeometric_exact_all_successes() {
assert_eq!(hyper_geometric_exactly(10, 10, 5, 5), 1.0); let input = &HyperGeometricInput::new(10, 10, 5, 5).unwrap();
assert_eq!(hyper_geometric_exactly(10, 10, 5, 4), 0.0); assert_eq!(hyper_geometric_exactly(input), 1.0);
let input = &HyperGeometricInput::new(10, 10, 5, 4).unwrap();
assert_eq!(hyper_geometric_exactly(input), 0.0);
} }
#[test] #[test]
fn test_hypergeometric_exact_no_successes() { fn test_hypergeometric_exact_no_successes() {
assert_eq!(hyper_geometric_exactly(10, 0, 5, 0), 1.0); let input = &HyperGeometricInput::new(10, 0, 5, 0).unwrap();
assert_eq!(hyper_geometric_exactly(10, 0, 5, 1), 0.0); assert_eq!(hyper_geometric_exactly(input), 1.0);
let input = &HyperGeometricInput::new(10, 0, 5, 1).unwrap();
assert_eq!(hyper_geometric_exactly(input), 0.0);
} }
#[test] #[test]
fn test_hypergeometric_exact() { fn test_hypergeometric_exact() {
assert_eq!(hyper_geometric_exactly(10, 3, 5, 2), 5.0 / 12.0); let input = &HyperGeometricInput::new(10, 3, 5, 2).unwrap();
assert_eq!(hyper_geometric_exactly(input), 5.0 / 12.0);
}
#[test]
fn test_hypergeometric_aces_poker() {
let input = HyperGeometricInput::new(52, 4, 5, 4).unwrap();
let exact = 1.846892603195124e-5;
assert_eq!(
hyper_geometric(input),
HyperGeometricProb {
exactly: exact,
less_than: 1.0 - exact,
less_or_equal: 1.0,
greater_than: 0.0,
greater_or_equal: exact
}
);
}
#[test]
fn test_binom_exact_all_success() {
let (p_powers, pc_powers) = powers(1.0, 5);
let input = &BinomialInput::new(1.0, 5, 5).unwrap();
assert_eq!(binomial_exactly(input, &p_powers, &pc_powers), 1.0);
let input = &BinomialInput::new(1.0, 5, 4).unwrap();
assert_eq!(binomial_exactly(input, &p_powers, &pc_powers), 0.0);
}
#[test]
fn test_binom_exact_no_success() {
let (p_powers, pc_powers) = powers(0.0, 5);
let input = &BinomialInput::new(0.0, 5, 0).unwrap();
assert_eq!(binomial_exactly(input, &p_powers, &pc_powers), 1.0);
let input = &BinomialInput::new(0.0, 5, 1).unwrap();
assert_eq!(binomial_exactly(input, &p_powers, &pc_powers), 0.0);
}
#[test]
fn test_binomial_exact() {
let (p_powers, pc_powers) = powers(0.5, 5);
let input = &BinomialInput::new(0.5, 5, 3).unwrap();
assert_eq!(binomial_exactly(input, &p_powers, &pc_powers), 10.0 / 32.0);
}
#[test]
fn test_binomial() {
// 10.0 / 32.0
let input = BinomialInput::new(0.5, 5, 3).unwrap();
assert_eq!(
binomial(input),
BinomialProb {
exactly: 10.0 / 32.0,
less_than: 16.0 / 32.0,
less_or_equal: 26.0 / 32.0,
greater_than: 6.0 / 32.0,
greater_or_equal: 16.0 / 32.0,
}
);
} }
} }

View file

@ -6,7 +6,7 @@ use leptos::prelude::{
use leptos::{IntoView, component, view}; use leptos::{IntoView, component, view};
use leptos_fluent::move_tr; use leptos_fluent::move_tr;
use crate::calc::{binomial, hyper_geometric}; use crate::calc::{BinomialInput, HyperGeometricInput, binomial, hyper_geometric};
#[component] #[component]
pub fn HyperCalculator() -> impl IntoView { pub fn HyperCalculator() -> impl IntoView {
@ -15,12 +15,13 @@ pub fn HyperCalculator() -> impl IntoView {
let (sample, set_sample) = signal(0u8); let (sample, set_sample) = signal(0u8);
let (sample_successes, set_sample_successes) = signal(0u8); let (sample_successes, set_sample_successes) = signal(0u8);
let result = move || { let result = move || {
hyper_geometric( HyperGeometricInput::new(
population.get(), population.get(),
successes.get(), successes.get(),
sample.get(), sample.get(),
sample_successes.get(), sample_successes.get(),
) )
.map(hyper_geometric)
.unwrap_or_default() .unwrap_or_default()
}; };
view! { view! {
@ -119,11 +120,12 @@ pub fn BinomCalculator() -> impl IntoView {
let (trials_number, set_trials_number) = signal(0u8); let (trials_number, set_trials_number) = signal(0u8);
let (successes_number, set_successes_number) = signal(0u8); let (successes_number, set_successes_number) = signal(0u8);
let result = move || { let result = move || {
binomial( BinomialInput::new(
success_probability.get(), success_probability.get(),
trials_number.get(), trials_number.get(),
successes_number.get(), successes_number.get(),
) )
.map(binomial)
.unwrap_or_default() .unwrap_or_default()
}; };
view! { view! {
@ -135,6 +137,7 @@ pub fn BinomCalculator() -> impl IntoView {
type="number" type="number"
min=0 min=0
max=1 max=1
prop:step=0.1
prop:value=success_probability prop:value=success_probability
on:input:target=move |ev| { on:input:target=move |ev| {
set_success_probability.set(ev.target().value().parse().unwrap_or_default()) set_success_probability.set(ev.target().value().parse().unwrap_or_default())