1212//! * Student's t
1313//! * Uniform
1414//! * Weighted Uniform
15+ //! * Log Normal
1516//! * There are two enums to represent probability distribution
1617//! * `OPDist<T>` : One parameter distribution (Bernoulli)
1718//! * `TPDist<T>` : Two parameter distribution (Uniform, Normal, Beta, Gamma)
226227//!
227228//! * Reference
228229//! * [Piecewise Rejection Sampling](https://axect.github.io/posts/006_prs/#22-weighted-uniform-distribution)
230+ //!
231+ //! ### Log Normal Distribution
232+ //!
233+ //! * Definition
234+ //! $$\text{LogNormal}(x | \mu, \sigma) = \frac{1}{x\sigma\sqrt{2\pi}} e^{-\frac{(\ln x -
235+ //! \mu)^2}{2\sigma^2}}$$
236+ //! where $\mu$ is the mean of the natural logarithm of the variable and $\sigma$ is the
237+ //! standard deviation of the natural logarithm of the variable.
238+ //! * Representative value
239+ //! * Mean: $e^{\mu + \frac{\sigma^2}{2}}$
240+ //! * Var: $(e^{\sigma^2} - 1)e^{2\mu + \sigma^2}$
241+ //! * To generate log-normal random samples, Peroxide uses the `rand_distr::LogNormal` distribution from the `rand_distr` crate.
229242
230243extern crate rand;
231244extern crate rand_distr;
@@ -267,6 +280,7 @@ pub enum TPDist<T: PartialOrd + SampleUniform + Copy + Into<f64>> {
267280 Normal ( T , T ) ,
268281 Beta ( T , T ) ,
269282 Gamma ( T , T ) ,
283+ LogNormal ( T , T ) ,
270284}
271285
272286pub struct WeightedUniform < T : PartialOrd + SampleUniform + Copy + Into < f64 > > {
@@ -312,7 +326,7 @@ impl WeightedUniform<f64> {
312326///
313327/// Ok(())
314328/// }
315- /// ```
329+ /// ```
316330pub fn new ( weights : Vec < f64 > , intervals : Vec < f64 > ) -> Result < Self > {
317331 let mut weights = weights;
318332 if weights. len ( ) == 0 {
@@ -497,6 +511,7 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> ParametricDist for TPDist
497511 Normal ( mu, sigma) => ( ( * mu) . into ( ) , ( * sigma) . into ( ) ) ,
498512 Beta ( a, b) => ( ( * a) . into ( ) , ( * b) . into ( ) ) ,
499513 Gamma ( a, b) => ( ( * a) . into ( ) , ( * b) . into ( ) ) ,
514+ LogNormal ( mu, sigma) => ( ( * mu) . into ( ) , ( * sigma) . into ( ) ) ,
500515 }
501516 }
502517}
@@ -683,7 +698,7 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for TPDist<T> {
683698 // }
684699 Gamma ( shape, scale) => {
685700 let gamma =
686- rand_distr:: Gamma :: < f64 > :: new ( ( * shape) . into ( ) , ( * scale) . into ( ) ) . unwrap ( ) ;
701+ rand_distr:: Gamma :: < f64 > :: new ( ( * shape) . into ( ) , 1f64 / ( * scale) . into ( ) ) . unwrap ( ) ;
687702 gamma. sample_iter ( rng) . take ( n) . collect ( )
688703 } // Gamma(a, b) => {
689704 // let a_f64 = (*a).into();
@@ -711,6 +726,11 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for TPDist<T> {
711726 // }
712727 // v
713728 // }
729+ LogNormal ( mu, sigma) => {
730+ let log_normal =
731+ rand_distr:: LogNormal :: < f64 > :: new ( ( * mu) . into ( ) , ( * sigma) . into ( ) ) . unwrap ( ) ;
732+ log_normal. sample_iter ( rng) . take ( n) . collect ( )
733+ }
714734 }
715735 }
716736
@@ -753,6 +773,16 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for TPDist<T> {
753773 * x. into ( ) . powf ( a_f64 - 1f64 )
754774 * E . powf ( -b_f64 * x. into ( ) )
755775 }
776+ LogNormal ( mu, sigma) => {
777+ let mu_f64 = ( * mu) . into ( ) ;
778+ let sigma_f64 = ( * sigma) . into ( ) ;
779+ if x. into ( ) <= 0f64 {
780+ 0f64
781+ } else {
782+ 1f64 / ( x. into ( ) * sigma_f64 * ( 2f64 * std:: f64:: consts:: PI ) . sqrt ( ) )
783+ * E . powf ( -( ( x. into ( ) . ln ( ) - mu_f64) . powi ( 2 ) ) / ( 2f64 * sigma_f64. powi ( 2 ) ) )
784+ }
785+ }
756786 }
757787 }
758788
@@ -791,6 +821,11 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> RNG for TPDist<T> {
791821
792822 inc_gamma ( a, b * x)
793823 }
824+ LogNormal ( mu, sigma) => {
825+ phi (
826+ ( x. ln ( ) - ( * mu) . into ( ) ) / ( * sigma) . into ( ) ,
827+ )
828+ }
794829 }
795830 }
796831}
@@ -886,6 +921,11 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> Statistics for TPDist<T>
886921 Normal ( m, _s) => ( * m) . into ( ) ,
887922 Beta ( a, b) => ( * a) . into ( ) / ( ( * a) . into ( ) + ( * b) . into ( ) ) ,
888923 Gamma ( a, b) => ( * a) . into ( ) / ( * b) . into ( ) ,
924+ LogNormal ( mu, sigma) => {
925+ let mu_f64 = ( * mu) . into ( ) ;
926+ let sigma_f64 = ( * sigma) . into ( ) ;
927+ E . powf ( mu_f64 + 0.5 * sigma_f64. powi ( 2 ) )
928+ }
889929 }
890930 }
891931
@@ -900,6 +940,11 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> Statistics for TPDist<T>
900940 a_f64 * b_f64 / ( ( a_f64 + b_f64) . powi ( 2 ) * ( a_f64 + b_f64 + 1f64 ) )
901941 }
902942 Gamma ( a, b) => ( * a) . into ( ) / ( * b) . into ( ) . powi ( 2 ) ,
943+ LogNormal ( mu, sigma) => {
944+ let mu_f64 = ( * mu) . into ( ) ;
945+ let sigma_f64 = ( * sigma) . into ( ) ;
946+ ( E . powf ( sigma_f64. powi ( 2 ) ) - 1f64 ) * E . powf ( 2f64 * mu_f64 + sigma_f64. powi ( 2 ) )
947+ }
903948 }
904949 }
905950
@@ -910,6 +955,7 @@ impl<T: PartialOrd + SampleUniform + Copy + Into<f64>> Statistics for TPDist<T>
910955 Normal ( _m, s) => ( * s) . into ( ) ,
911956 Beta ( _a, _b) => self . var ( ) . sqrt ( ) ,
912957 Gamma ( _a, _b) => self . var ( ) . sqrt ( ) ,
958+ LogNormal ( _mu, _sigma) => self . var ( ) . sqrt ( ) ,
913959 }
914960 }
915961
0 commit comments