Axect / Peroxide

Rust numeric library with R, MATLAB & Python syntax
https://crates.io/crates/peroxide
Apache License 2.0
509 stars 31 forks source link

Optimization with uncertainties #26

Closed hombit closed 4 years ago

hombit commented 4 years ago

Is there a way to use Optimizer with weights? It would be very useful if your data include uncertainties sigma, a user can easily found weights as 1 / sigma^2.

If there is no such way, would you accept a PR where loss function is sum (f_i - y(x_i))^2 / sigma_i^2, Optimizer has additional field weight: Optional<Vec<f64>> and Optimizer::new accepts two or three column data (additional optional column is weight)?

Axect commented 4 years ago

Thanks to your interests. Although there is no official support for weight, but you can capture weight vector via closure.

For example, let's do linear regression with additional error term.

#[macro_use]
extern crate peroxide;
use peroxide::fuga::*;

fn main() {
    // Generate Data
    let x = seq(1, 10, 1);
    let eps_true = rnorm!(10, 0, 1);
    let y = x.fmap(|t| 2f64 * t + 1f64).add_v(&eps_true);
    let data = cbind(x.into(), y.into());

    // Generate Normal error (eps ~ N(0, 1^2))
    let eps_guess = rnorm!(10, 0, 1);

    // Optimize with new error
    let mut opt = Optimizer::new(data, |x, w| linear_regression(x, w, &eps_guess));
    let p = opt.set_init_param(c!(1, 1))
        .set_max_iter(50)
        .set_method(LevenbergMarquardt)
        .optimize();
    p.print();
    opt.get_error().print();
}

// Linear regression with error: y_i = w_0 + w_1 * x_i + epsilon_i
fn linear_regression(x: &Vec<f64>, w: Vec<Number>, epsilon: &Vec<f64>) -> Option<Vec<Number>> {
    Some(
        x.iter()
            .map(|&t| Number::from_f64(t))
            .zip(epsilon.iter().map(|&t| Number::from_f64(t)))
            .map(|(t, eps)| w[0] + w[1] * t + eps)
            .collect()
    )
}

Maybe add additional argument can make some confusion. For this reason, it is cautious to accept your good suggestion.

If you have more question or suggestion, then please let me know.