More general Exp

[saona-raimundo]

Background

What is your motivation?
I am implementing stochastic processes, continuous markov chains in particular, and I need an Exponential clock. I wanted to use rand_distr::Exp as an internal field, but I have the following two problems:

1. To include absorbent states, I need a zero rate exponential, i.e. that sample always INFINITY.
2. It would be nice to have a method set_lambda to change it, instead of having to construct a new one.

What type of application is this?
Numerical simulation

Feature request

Generalize the construction of Exp to include the case of lambda = 0.0.
Also, implement a "setter" method to change the lambda in an already created Exp.

I know that allowing `lambda = 0.0` could imply a comparison before each sample (checking if lambda is finite or not) and I do not know of any solution for this. If there is no way around to maintain performance, then I would prefer this feature not to be implemented.

[vks]

I'm not sure I understand what the problem is with constructing a new one. Is this about performance or convenience?

If you need this edge case behavior, couldn't you implement your own distribution on top of Exp?

[dhardy]

1. I don't understand whether this infinite case is merely a possible random sample (in which case it should happen rarely enough it probably doesn't matter) or is explicitly selected. But I understand your point: from this point of view, λ = 0 is a valid parametrisation. Perhaps we should add a second constructor fn new_allowing_zero(lambda: N) -> Self?
2. Using distr = Exp::new(...) should have the same effect (probably exactly the same code when optimised). This is one reason we always mark constructors with #[inline].

[saona-raimundo]

Thank you for the quick answer!!

Sorry for not being clear. I will try again.

General parametrization

λ = 0 is a valid parametrisation.

Yes, from my point of view it is (Wikipedia does not think the same).

Implications

Perhaps we should add a second constructor fn new_allowing_zero(lambda: N) -> Self?

The problem is that the implementation changes. An exponential with λ = 0 is always (samples always) infinity. But, an exponential with λ > 0 samples always finite and has the same distribution as Exp(1) / λ.
Therefore, I thought that allowing λ = 0 in the construction and using the naive implementation (use of if for every sample) will be "less performant for the general use" since the variable Exp is usually used with λ > 0.
If that is the case and there is no way around it,

couldn't you implement your own distribution on top of Exp?

I am okay using my own wrapper, I just wanted to ask if you had better ideas.

Setter method

I'm not sure I understand what the problem is with constructing a new one.

There is no problem with constructing a new one, a method set_lambda would be purely for convenience. Feel free to say no, I just thought the code would look nicer.

[dhardy]

(use of if for every sample)

I was thinking optimise for the common case, waste a sample and multiply by INF in the zero case. Is that not acceptable?

purely for convenience

When designing an API for general usage, it is always a good idea to minimise complexity: avoid methods you don't need. And I don't think a setter adds enough.

[saona-raimundo]

waste a sample and multiply by INF in the zero case. Is that not acceptable?

Good!! I didn't think about it. In fact, there is no need to waste a sample! This is the current implementation.

#[derive(Clone, Copy, Debug)]
pub struct Exp<N> {
    /// `lambda` stored as `1/lambda`, since this is what we scale by.
    lambda_inverse: N
}

/// Error type returned from `Exp::new`.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Error {
    /// `lambda <= 0` or `nan`.
    LambdaTooSmall,
}

impl<N: Float> Exp<N>
where Exp1: Distribution<N>
{
    /// Construct a new `Exp` with the given shape parameter
    /// `lambda`.
    #[inline]
    pub fn new(lambda: N) -> Result<Exp<N>, Error> {
        if !(lambda > N::from(0.0)) {
            return Err(Error::LambdaTooSmall);
        }
        Ok(Exp { lambda_inverse: N::from(1.0) / lambda })
    }
}

impl<N: Float> Distribution<N> for Exp<N>
where Exp1: Distribution<N>
{
    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> N {
        rng.sample(Exp1) * self.lambda_inverse
    }
}

One could change the constructor method (or create another constructor to not have breaking change) like this.

    pub fn new(lambda: N) -> Result<Exp<N>, Error> {
        if lambda > N::from(0.0) {
            return Ok(Exp { lambda_inverse: N::from(1.0) / lambda });
        } else if lambda < N::from(0.0) {
            return Err(Error::LambdaTooSmall);
        } else {
        Ok(Exp { lambda_inverse: N::infinity() })
        }
    }
}

Consequences

1. No performance change for the common use case.
2. Ask to implement infinity for the rand_distr::Float trait.
3. Change documentation of LambdaTooSmall for strictly negative lambda.

[dhardy]

@rasa200 that's roughly what I was thinking, yes. It does still sample Exp1 when calling sample(), but I think that's acceptable enough.

@vks any thoughts on whether we should allow zero-rate in the new constructor? The only issue is whether some users expect an error in this case. It may be okay just to relax the bound to lambda >= 0?

[vks]

I think either way is fine as long as it is documented.

[saona-raimundo]

It does still sample Exp1 when calling sample(), but I think that's acceptable enough.

Sure! It is good!

[dhardy]

Okay then, I guess the only code change we need is to Exp::new. Documentation-wise, we need a clear note on that constructor, tweaks to the error type, and a note in the changelog.

@rasa200 would you like to make a PR?

[saona-raimundo]

@dhardy Happy to make a PR :D It will be my first ever.

One question: To use N::infinity(), I would need to add the infinity method to the rand_distr::Float trait. Is that okay?

Another possibility is using
N::from(f64::INFINITY) instead if N::infinity(), which:

• Does not require to change the Float trait
• Ask for possibly one unnecessary casting in the "non-common" case.

Historical note: In rand, version 0.7.3, rand::distributions::Exp can be created only with f64

[dhardy]

I would need to add the infinity method to the rand_distr::Float trait. Is that okay?

You shouldn't need to. I just tested:

    /// Construct a new `Exp` with the given shape parameter
    /// `lambda`.
    #[inline]
    pub fn new(lambda: N) -> Result<Exp<N>, Error> {
        if !(lambda >= N::from(0.0)) {
            return Err(Error::LambdaTooSmall);
        }
        Ok(Exp {
            lambda_inverse: N::from(1.0) / lambda,
        })
    }
}

...

#[cfg(test)]
mod test {
    use super::*;

    #[test]
    fn test_zero() {
        let d = Exp::new(0.0).unwrap();
        assert_eq!(d.sample(&mut crate::test::rng(21)), std::f64::INFINITY);
    }

---

Yes, supporting both float types is a recent change.

[saona-raimundo]

True!! My bad, I should have tried first.

[saona-raimundo]

There is a PR :D

I had a bit of travel with the documentation (see PR), because of the possibility of a custom implementation of 1 / 0 for types implementing the rand_distr::Float trait. Other than that, I should have taken all the discussion into account.