emcee

在 Rust 中重新实现 Dan Foreman-Mackey 的 emcee。

请参阅此处的主文档

fitting_a_model_to_data 示例 是从 emcee 文档中 "fitting a model to data" 示例的重新创建。

归属

如果你在工作中使用了 emcee，请引用 Dan 的论文（《arXiv》, 《ADS》, 《BibTeX》）。

原始 MIT 许可证的副本在 DFM-LICENSE 中。

基本用法

实现模型

采样器需要一个实现了 emcee::Prob 的结构体，例如

use emcee::{Guess, Prob};

struct Model;

impl Prob for Model {
    fn lnlike(&self, params: &Guess) -> f32 {
        // Insert actual implementation here
        0f32
    }

    fn lnprior(&self, params: &Guess) -> f32 {
        // Insert actual implementation here
        0f32
    }
}

该特性为 lnprob 提供了一个默认实现，该实现根据贝叶斯定理计算似然函数和先验概率的乘积（对数空间中的和）。无效的先验值通过从先验函数中返回 -std::f32::INFINITY 来标记。请注意，您的实现可能需要外部数据。这些数据应包含在您的 Model 类中，例如

struct Model<'a> {
    x: &'a [f32],
    y: &'a [f32],
}

// Linear model y = m * x + c
impl<'a> Prob for Model<'a> {
    fn lnlike(&self, params: &Guess) -> f32 {
        let m = params[0];
        let c = params[1];

        -0.5 * self.x.iter().zip(self.y)
            .map(|(xval, yval)| {
                let model = m * xval + c;
                let residual = (yval - model).powf(2.0);
                residual
            }).sum::<f32>()
    }

    fn lnprior(&self, params: &Guess) -> f32 {
        // unimformative priors
        0.0f32
    }
}

初始猜测

接下来，构建一个初始猜测。一个 Guess 表示一个建议参数向量

use emcee::Guess;

let initial_guess = Guess::new(&[0.0f32, 0.0f32]);

此create实现的采样器使用多个walker，因此每个walker都需要复制一次初始猜测，通常从初始位置散布开来，以帮助探索问题参数空间。这可以通过create_initial_guess方法实现。

let nwalkers = 100;
let perturbed_guess = initial_guess.create_initial_guess(nwalkers);
assert_eq!(perturbed_guess.len(), nwalkers);

构建采样器

采样器生成新的参数向量，使用用户提供的概率模型评估概率，接受更可能的参数向量，并进行多次迭代。

采样器需要知道要使用的walker数量，这必须是一个偶数，并且至少是参数向量大小的两倍。它还需要你的参数向量的大小，以及你的概率结构（实现了Prob）

let nwalkers = 100;
let ndim = 2;  // m and c

// Build a linear model y = m * x + c (see above)

let initial_x = [0.0f32, 1.0f32, 2.0f32];
let initial_y = [5.0f32, 7.0f32, 9.0f32];

let model = Model {
    x: &initial_x,
    y: &initial_y,
};

let sampler = emcee::EnsembleSampler::new(nwalkers, ndim, &model)
    .expect("could not create sampler");

然后运行采样器

let niterations = 100;
sampler.run_mcmc(&perturbed_guess, niterations).expect("error running sampler");

迭代采样

有时在采样器的每个提议步骤中获取提议和评估的内部值是有用的。在Python版本中，sample方法是一个生成器，可以迭代以评估样本步骤。

在这个Rust版本中，我们通过公开sample方法来提供此功能，该方法接受一个回调，每次迭代都会调用一次，传递一个单一的Step对象。例如

sampler.sample(&perturbed_guess, niterations, |step| {
    println!("Current iteration: {}", step.iteration);
    println!("Current guess vectors: {:?}", step.pos);
    println!("Current log posterior probabilities: {:?}", step.lnprob);
});

研究结果

样本存储在采样器的flatchain中，这是通过采样器上的flatchain方法构建的

let flatchain = sampler.flatchain().unwrap();

for (i, guess) in flatchain.iter().enumerate() {
    // Skip possible "burn-in" phase
    if i < 50 * nwalkers {
        continue;
    }

    println!("Iteration {}; m={}, c={}", i, guess[0], guess[1]);
}