Rust version of https://gist.github.com/TanyaS08/5ac14d4969bd5d4b772bcdc533999c8f

sr2ch3h1.rs

/// You need to add these lines to Cargo.toml's [dependencies]:
/// statrs = "0.15.0"
/// nalgebra = "0.29.0"

use statrs::distribution::{Binomial, Discrete};
use nalgebra::base::DVector;
use nalgebra::dvector;

// Builds a vector of evenly spaced floats within a closed interval [a, b].
pub fn linspace(a: f64, b: f64, n: usize) -> DVector<f64> {
    let step_size = (b - a) / (n - 1) as f64;
    DVector::from_fn(n, |i, _| a + (i as f64) * step_size)
}

fn main() {
    let birth1 = dvector![1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0,0,0,0,1,1,1,0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,0,1,1,1,1,1,0,0,1,0,1,1,0,1,0,1,1,1,0,1,1,1,1];
    let birth2 = dvector![0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0,1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,1,0,0,0,1,1,1,0,0,0,0];

    let grid_size = 8128; // Perfect!
    let ps = linspace(0.0, 1.0, grid_size);

    let prior = DVector::repeat(grid_size, 1.0);

    let boys = birth1.sum() + birth2.sum();
    let n = (birth1.len() + birth2.len()) as u64;
    //let girls = n - boys;

    let likelihood = ps.map(|p|
        Binomial::new(p, n).unwrap().pmf(boys)
    );

    //let upost = likelihood.component_mul(&prior);
    let upost = likelihood.zip_map(&prior,
        |x: f64, y: f64| (x.ln() + y.ln()).exp()
    );

    println!("{}", ps[upost.imax()]); // imax() returns the index of the maximum.
}