diff --git a/examples/3d-plot2.rs b/examples/3d-plot2.rs index f2704658..38800759 100644 --- a/examples/3d-plot2.rs +++ b/examples/3d-plot2.rs @@ -1,11 +1,20 @@ use plotters::prelude::*; +use std::f64::consts::PI; fn pdf(x: f64, y: f64) -> f64 { - const SDX: f64 = 0.1; - const SDY: f64 = 0.1; - const A: f64 = 5.0; - let x = x as f64 / 10.0; - let y = y as f64 / 10.0; - A * (-x * x / 2.0 / SDX / SDX - y * y / 2.0 / SDY / SDY).exp() + // see https://mathworld.wolfram.com/BivariateNormalDistribution.html + // assumes x_bar = 0 and y_bar = 0 + const SDX: f64 = 0.3; + const SDY: f64 = 0.3; + const vx: f64 = SDX * SDX; + const vy : f64 = SDY * SDY; + const sdx_sdy: f64 = SDX * SDY; + const RHO : f64 = 0.9; + const c1: f64 = 1.0 - RHO * RHO; + let c2 : f64 = c1.sqrt(); + let A: f64 = 1.0 / (2.0 * PI * sdx_sdy * c2); + let Z: f64 = x * x / vx - 2.0 * RHO * x * y / sdx_sdy + y * y / vy; + let res = A * (-Z / 2.0 / c1).exp(); + res } fn main() -> Result<(), Box> { @@ -17,7 +26,7 @@ fn main() -> Result<(), Box> { let mut chart = ChartBuilder::on(&root) .caption("2D Guassian PDF", ("sans-serif", 20)) - .build_cartesian_3d(-3.0..3.0, 0.0..6.0, -3.0..3.0)?; + .build_cartesian_3d(-1.0..1.0, 0.0..3.0, -1.0..1.0)?; chart.with_projection(|mut p| { p.pitch = 1.57 - (1.57 - pitch as f64 / 50.0).abs(); p.scale = 0.7; @@ -28,8 +37,8 @@ fn main() -> Result<(), Box> { chart.draw_series( SurfaceSeries::xoz( - (-15..=15).map(|x| x as f64 / 5.0), - (-15..=15).map(|x| x as f64 / 5.0), + (-30..=30).map(|x| x as f64 / 30.0), + (-30..=30).map(|x| x as f64 / 30.0), pdf, ) .style_func(&|&v| {