|
| 1 | +// Copyright 2021 Developers of the Rand project. |
| 2 | +// |
| 3 | +// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or |
| 4 | +// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license |
| 5 | +// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your |
| 6 | +// option. This file may not be copied, modified, or distributed |
| 7 | +// except according to those terms. |
| 8 | + |
| 9 | +//! The Skew Normal distribution. |
| 10 | +
|
| 11 | +use crate::{Distribution, StandardNormal}; |
| 12 | +use core::fmt; |
| 13 | +use num_traits::Float; |
| 14 | +use rand::Rng; |
| 15 | + |
| 16 | +/// The [skew normal distribution] `SN(location, scale, shape)`. |
| 17 | +/// |
| 18 | +/// The skew normal distribution is a generalization of the |
| 19 | +/// [`Normal`] distribution to allow for non-zero skewness. |
| 20 | +/// |
| 21 | +/// It has the density function, for `scale > 0`, |
| 22 | +/// `f(x) = 2 / scale * phi((x - location) / scale) * Phi(alpha * (x - location) / scale)` |
| 23 | +/// where `phi` and `Phi` are the density and distribution of a standard normal variable. |
| 24 | +/// |
| 25 | +/// # Example |
| 26 | +/// |
| 27 | +/// ``` |
| 28 | +/// use rand_distr::{SkewNormal, Distribution}; |
| 29 | +/// |
| 30 | +/// // location 2, scale 3, shape 1 |
| 31 | +/// let skew_normal = SkewNormal::new(2.0, 3.0, 1.0).unwrap(); |
| 32 | +/// let v = skew_normal.sample(&mut rand::thread_rng()); |
| 33 | +/// println!("{} is from a SN(2, 3, 1) distribution", v) |
| 34 | +/// ``` |
| 35 | +/// |
| 36 | +/// # Implementation details |
| 37 | +/// |
| 38 | +/// We are using the algorithm from [A Method to Simulate the Skew Normal Distribution]. |
| 39 | +/// |
| 40 | +/// [skew normal distribution]: https://en.wikipedia.org/wiki/Skew_normal_distribution |
| 41 | +/// [`Normal`]: struct.Normal.html |
| 42 | +/// [A Method to Simulate the Skew Normal Distribution]: |
| 43 | +/// Ghorbanzadeh, D. , Jaupi, L. and Durand, P. (2014) |
| 44 | +/// [A Method to Simulate the Skew Normal Distribution](https://dx.doi.org/10.4236/am.2014.513201). |
| 45 | +/// Applied Mathematics, 5, 2073-2076. |
| 46 | +#[derive(Clone, Copy, Debug)] |
| 47 | +#[cfg_attr(feature = "serde1", derive(serde::Serialize, serde::Deserialize))] |
| 48 | +pub struct SkewNormal<F> |
| 49 | +where |
| 50 | + F: Float, |
| 51 | + StandardNormal: Distribution<F>, |
| 52 | +{ |
| 53 | + location: F, |
| 54 | + scale: F, |
| 55 | + shape: F, |
| 56 | +} |
| 57 | + |
| 58 | +/// Error type returned from `SkewNormal::new`. |
| 59 | +#[derive(Clone, Copy, Debug, PartialEq, Eq)] |
| 60 | +pub enum Error { |
| 61 | + /// The scale parameter is not finite or it is less or equal to zero. |
| 62 | + ScaleTooSmall, |
| 63 | + /// The shape parameter is not finite. |
| 64 | + BadShape, |
| 65 | +} |
| 66 | + |
| 67 | +impl fmt::Display for Error { |
| 68 | + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { |
| 69 | + f.write_str(match self { |
| 70 | + Error::ScaleTooSmall => { |
| 71 | + "scale parameter is either non-finite or it is less or equal to zero in skew normal distribution" |
| 72 | + } |
| 73 | + Error::BadShape => "shape parameter is non-finite in skew normal distribution", |
| 74 | + }) |
| 75 | + } |
| 76 | +} |
| 77 | + |
| 78 | +#[cfg(feature = "std")] |
| 79 | +#[cfg_attr(doc_cfg, doc(cfg(feature = "std")))] |
| 80 | +impl std::error::Error for Error {} |
| 81 | + |
| 82 | +impl<F> SkewNormal<F> |
| 83 | +where |
| 84 | + F: Float, |
| 85 | + StandardNormal: Distribution<F>, |
| 86 | +{ |
| 87 | + /// Construct, from location, scale and shape. |
| 88 | + /// |
| 89 | + /// Parameters: |
| 90 | + /// |
| 91 | + /// - location (unrestricted) |
| 92 | + /// - scale (must be finite and larger than zero) |
| 93 | + /// - shape (must be finite) |
| 94 | + #[inline] |
| 95 | + pub fn new(location: F, scale: F, shape: F) -> Result<SkewNormal<F>, Error> { |
| 96 | + if !scale.is_finite() || !(scale > F::zero()) { |
| 97 | + return Err(Error::ScaleTooSmall); |
| 98 | + } |
| 99 | + if !shape.is_finite() { |
| 100 | + return Err(Error::BadShape); |
| 101 | + } |
| 102 | + Ok(SkewNormal { |
| 103 | + location, |
| 104 | + scale, |
| 105 | + shape, |
| 106 | + }) |
| 107 | + } |
| 108 | + |
| 109 | + /// Returns the location of the distribution. |
| 110 | + pub fn location(&self) -> F { |
| 111 | + self.location |
| 112 | + } |
| 113 | + |
| 114 | + /// Returns the scale of the distribution. |
| 115 | + pub fn scale(&self) -> F { |
| 116 | + self.scale |
| 117 | + } |
| 118 | + |
| 119 | + /// Returns the shape of the distribution. |
| 120 | + pub fn shape(&self) -> F { |
| 121 | + self.shape |
| 122 | + } |
| 123 | +} |
| 124 | + |
| 125 | +impl<F> Distribution<F> for SkewNormal<F> |
| 126 | +where |
| 127 | + F: Float, |
| 128 | + StandardNormal: Distribution<F>, |
| 129 | +{ |
| 130 | + fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F { |
| 131 | + let linear_map = |x: F| -> F { x * self.scale + self.location }; |
| 132 | + let u_1: F = rng.sample(StandardNormal); |
| 133 | + if self.shape == F::zero() { |
| 134 | + linear_map(u_1) |
| 135 | + } else { |
| 136 | + let u_2 = rng.sample(StandardNormal); |
| 137 | + let (u, v) = (u_1.max(u_2), u_1.min(u_2)); |
| 138 | + if self.shape == -F::one() { |
| 139 | + linear_map(v) |
| 140 | + } else if self.shape == F::one() { |
| 141 | + linear_map(u) |
| 142 | + } else { |
| 143 | + let normalized = ((F::one() + self.shape) * u + (F::one() - self.shape) * v) |
| 144 | + / ((F::one() + self.shape * self.shape).sqrt() |
| 145 | + * F::from(core::f64::consts::SQRT_2).unwrap()); |
| 146 | + linear_map(normalized) |
| 147 | + } |
| 148 | + } |
| 149 | + } |
| 150 | +} |
| 151 | + |
| 152 | +#[cfg(test)] |
| 153 | +mod tests { |
| 154 | + use super::*; |
| 155 | + |
| 156 | + fn test_samples<F: Float + core::fmt::Debug, D: Distribution<F>>( |
| 157 | + distr: D, zero: F, expected: &[F], |
| 158 | + ) { |
| 159 | + let mut rng = crate::test::rng(213); |
| 160 | + let mut buf = [zero; 4]; |
| 161 | + for x in &mut buf { |
| 162 | + *x = rng.sample(&distr); |
| 163 | + } |
| 164 | + assert_eq!(buf, expected); |
| 165 | + } |
| 166 | + |
| 167 | + #[test] |
| 168 | + #[should_panic] |
| 169 | + fn invalid_scale_nan() { |
| 170 | + SkewNormal::new(0.0, core::f64::NAN, 0.0).unwrap(); |
| 171 | + } |
| 172 | + |
| 173 | + #[test] |
| 174 | + #[should_panic] |
| 175 | + fn invalid_scale_zero() { |
| 176 | + SkewNormal::new(0.0, 0.0, 0.0).unwrap(); |
| 177 | + } |
| 178 | + |
| 179 | + #[test] |
| 180 | + #[should_panic] |
| 181 | + fn invalid_scale_negative() { |
| 182 | + SkewNormal::new(0.0, -1.0, 0.0).unwrap(); |
| 183 | + } |
| 184 | + |
| 185 | + #[test] |
| 186 | + #[should_panic] |
| 187 | + fn invalid_scale_infinite() { |
| 188 | + SkewNormal::new(0.0, core::f64::INFINITY, 0.0).unwrap(); |
| 189 | + } |
| 190 | + |
| 191 | + #[test] |
| 192 | + #[should_panic] |
| 193 | + fn invalid_shape_nan() { |
| 194 | + SkewNormal::new(0.0, 1.0, core::f64::NAN).unwrap(); |
| 195 | + } |
| 196 | + |
| 197 | + #[test] |
| 198 | + #[should_panic] |
| 199 | + fn invalid_shape_infinite() { |
| 200 | + SkewNormal::new(0.0, 1.0, core::f64::INFINITY).unwrap(); |
| 201 | + } |
| 202 | + |
| 203 | + #[test] |
| 204 | + fn valid_location_nan() { |
| 205 | + SkewNormal::new(core::f64::NAN, 1.0, 0.0).unwrap(); |
| 206 | + } |
| 207 | + |
| 208 | + #[test] |
| 209 | + fn skew_normal_value_stability() { |
| 210 | + test_samples( |
| 211 | + SkewNormal::new(0.0, 1.0, 0.0).unwrap(), |
| 212 | + 0f32, |
| 213 | + &[-0.11844189, 0.781378, 0.06563994, -1.1932899], |
| 214 | + ); |
| 215 | + test_samples( |
| 216 | + SkewNormal::new(0.0, 1.0, 0.0).unwrap(), |
| 217 | + 0f64, |
| 218 | + &[ |
| 219 | + -0.11844188827977231, |
| 220 | + 0.7813779637772346, |
| 221 | + 0.06563993969580051, |
| 222 | + -1.1932899004186373, |
| 223 | + ], |
| 224 | + ); |
| 225 | + test_samples( |
| 226 | + SkewNormal::new(core::f64::INFINITY, 1.0, 0.0).unwrap(), |
| 227 | + 0f64, |
| 228 | + &[ |
| 229 | + core::f64::INFINITY, |
| 230 | + core::f64::INFINITY, |
| 231 | + core::f64::INFINITY, |
| 232 | + core::f64::INFINITY, |
| 233 | + ], |
| 234 | + ); |
| 235 | + test_samples( |
| 236 | + SkewNormal::new(core::f64::NEG_INFINITY, 1.0, 0.0).unwrap(), |
| 237 | + 0f64, |
| 238 | + &[ |
| 239 | + core::f64::NEG_INFINITY, |
| 240 | + core::f64::NEG_INFINITY, |
| 241 | + core::f64::NEG_INFINITY, |
| 242 | + core::f64::NEG_INFINITY, |
| 243 | + ], |
| 244 | + ); |
| 245 | + } |
| 246 | + |
| 247 | + #[test] |
| 248 | + fn skew_normal_value_location_nan() { |
| 249 | + let skew_normal = SkewNormal::new(core::f64::NAN, 1.0, 0.0).unwrap(); |
| 250 | + let mut rng = crate::test::rng(213); |
| 251 | + let mut buf = [0.0; 4]; |
| 252 | + for x in &mut buf { |
| 253 | + *x = rng.sample(&skew_normal); |
| 254 | + } |
| 255 | + for value in buf.iter() { |
| 256 | + assert!(value.is_nan()); |
| 257 | + } |
| 258 | + } |
| 259 | +} |
0 commit comments