Generalized inverse Gaussian

patwa67 commented 7 years ago

Someone that plans to take on the task of implementing the generalized inverse Gaussian distribution? DOI: 10.1007/s11222-013-9387-3 https://cran.r-project.org/web/packages/GIGrvg/GIGrvg.pdf

ararslan commented 7 years ago

I don't know that anyone here is already planning to take that on, so a PR would be welcome.

However, do note that the R code to which you linked it licensed under GPL. That means that you can't simply translate the code to Julia, since that would constitute a "derived work," and would require the Julia code to be licensed under GPL as well. We wouldn't be able to accept the PR, since this package is licensed under MIT.

Perhaps you already know all that, just figured I'd mention it just in case. 😉

patwa67 commented 7 years ago

I don't think I'm the right person to lead a PR for this as Julia is very new to me. I would be happy to follow the process if someone else could take the lead. There is very useful pseudo-code in the paper, but it would take weeks for me to get that into Julia.

ararslan commented 7 years ago

If it's something you'd be interested in doing then I encourage you to try your hand at it. I find things like this can be incredibly rewarding to work on, and if you need any guidance along the way I'd be happy to lend a hand. Up to you though; if you'd rather leave it to someone else that's of course totally fine as well. I just don't want you to feel like you can't or shouldn't just because you're new. Everyone starts somewhere!

patwa67 commented 7 years ago

I can give it a try, but would need help on how to start up the PR. Could you shortly outline how to do that?

ararslan commented 7 years ago

Of course! Here's my personally preferred workflow:

Fork the repository to your account (click Fork in the upper right corner)
Go somewhere on your computer, git clone https://github.com/JuliaStats/Distributions.jl.git
git remote add patwa67 https://github.com/patwa67/Distributions.jl.git (this is your fork)
git checkout -b yourbranch (where yourbranch is whatever you want to call it)
Create a file src/univariate/continuous/generalizedinversegaussian.jl
Create the distribution's type and define methods
- Here we define a distribution as an immutable that inherits from an abstract type.
- See src/univariate/continuous/inversegaussian.jl for an example of the interface required for a distribution (things like pdf, cdf, rand, etc.)
- Make sure to add doc strings (again see other distributions for examples)
git add src/univariate/continuous/generalizedinversegaussian.jl (stage the file for commit)
Create a file test/generalizedinversegaussian.jl and test all of the functionality you've implemented for your distribution
- See normalinversegaussian.jl for an example
git add test/generalizedinversegaussian.jl
Add documentation in the appropriate ReST file(s) in doc/source
git add doc/source/(any files you may have created)
git commit -m "Add generalized inverse gaussian"
git push patwa67 yourbranch
Go to the main page for this repo. There will be a yellow bar from which you can select "open PR."
We can go from there. 🙂

That list may look daunting, but I promise it's less complicated than it looks. @andreasnoack, have I omitted anything here? Is there anything else patwa should know to get started?

patwa67 commented 7 years ago

I did start on this issue, but realized that I won't have time to finish this in time. I have a start template for generalizedinversegaussian.jl and would be happy if someone else could take over this task and open a proper PR. Just let me know who I should send the file to.

ararslan commented 7 years ago

finish this in time

There's no rush. 😉 But if you don't want to continue work on it, you could just have a commit for your preliminary changes on your fork, then whoever wants to pick up where you left off can use the commit from your fork.

patwa67 commented 7 years ago

I'm completely lost in this process. The only thing I have is a generalizedinversegaussian.jl (somewhere), but I doubt that it is on Github. I think I created a branch under master called GIG, but I cannot find it. It would be best if I just could send the file to someone.

ararslan commented 7 years ago

Git can be quite frustrating at times. However, the advantage of doing it through Git is that authorship information is retained; if done properly by the PR author, the PR will include your commit with you still shown as the commit author. If you're willing to give it another shot, try this:

You have a fork of Distributions on your account. That's great. That's step 1.
Navigate to the folder where you have Distributions cloned locally.
Ensure your fork is listed as a remote by submitting git remote -v. I'll assume now that the remote is registered as patwa67. If it's not in there, see the remote add step in my previous instructions.
Check out your generalized inverse Gaussian branch: git checkout GIG.
Ensure your changes are committed using git status (git add anything that isn't staged, git commit anything that isn't committed) and git log --oneline -n 10 (do the commits show up in the history?).
Push to your remote: git push patwa67 GIG:GIG.
Create a pull request here. Make sure to check "allow edits from maintainers."

That's it! Whoever takes it from here can work directly from your preliminary PR, plus you're able to contribute more to it as well.

Does that seem doable? Don't hestitate to let me know if you run into any problems or have questions.

patwa67 commented 7 years ago

I'm afraid that it doesn't work. The only thing a get is: patwa67 merged 1 commit into master from patwa67-GIG 3 minutes ago

but that doesn't turn up as a pull request.

ararslan commented 7 years ago

Does the commit exist on your GIG branch?

patwa67 commented 7 years ago

GitHub is not for me, here's the code:

doc"""
    GeneralizedInverseGaussian(a,b,p)

The *generalized inverse Gaussian distribution* with parameters a > 0, b > 0, p real,
and modified Bessel function of the second kind K_p, has probability density function

$f(x; a, b, p) = \frac{(a/b)^{p/2}}{2K_p(\sqrt{ab})}x^{(p-1)}
e^{-(ax+b/x)/2}, \quad x > 0$

``julia
GeneralizedInverseGaussian(a, b, p)    # Generalized Inverse Gaussian distribution with parameters parameters a > 0, b > 0 and p real

params(d)           # Get the parameters, i.e. (a, b, p)
``

External links

* [Generalized Inverse Gaussian distribution on Wikipedia](https://en.wikipedia.org/wiki/Generalized_inverse_Gaussian_distribution)

"""
immutable GeneralizedInverseGaussian{T<:Real} <: ContinuousUnivariateDistribution
    a::T
    b::T
    p::T

    function GeneralizedInverseGaussian(a::T, b::T, p::T)
        @check_args(GeneralizedInverseGaussian, a > zero(a) && b > zero(b))
        new(a, b, p)
    end
end

GeneralizedInverseGaussian{T<:Real}(a::T, b::T, p::T) = GeneralizedInverseGaussian{T}(a, b, p)
GeneralizedInverseGaussian(a::Real, b::Real, p::Real) = GeneralizedInverseGaussian(promote(a, b, p)...)
GeneralizedInverseGaussian(a::Integer, b::Integer, p::Integer) = GeneralizedInverseGaussian(Float64(a), Float64(b), Float64(p))

@distr_support GeneralizedInverseGaussian 0.0 Inf

#### Conversions

function convert{T <: Real, S <: Real}(::Type{GeneralizedInverseGaussian{T}}, a::S, b::S, p::S)
    GeneralizedInverseGaussian(T(a), T(b), T(p))
end
function convert{T <: Real, S <: Real}(::Type{GeneralizedInverseGaussian{T}}, d::GeneralizedInverseGaussian{S})
    GeneralizedInverseGaussian(T(d.a), T(d.b), T(d.p))
end

#### Parameters

params(d::GeneralizedInverseGaussian) = (d.a, d.b, d.p)
@inline partype{T<:Real}(d::InverseGaussian{T}) = T

#### Statistics

mean(d::GeneralizedInverseGaussian) = (sqrt(d.b) * besselk((d.p + 1), sqrt(d.a * d.b))) / (sqrt(d.a) * besselk(d.p, sqrt(d.a * d.b)))

var(d::GeneralizedInverseGaussian) = (b / a) * ((besselk(d.p + 2, sqrt(d.a * d.b)) / (besselk(d.p, sqrt(d.a * d.b)))) -
 (besselk(d.p + 1, sqrt(d.a * d.b)) / (besselk(d.p, sqrt(d.a * d.b))))^2)

mode(d::GeneralizedInverseGaussian) = ((d.p - 1) + sqrt((d.p - 1)^2 + d.a * d.b)) / d.a

#### Evaluation

function pdf{T<:Real}(d::GeneralizedInverseGaussian{T}, x::Real)
    if x > 0
        a, b, p = params(d)
        return (((a / b)^(p / 2)) / (2 * besselk(p, sqrt(a * b)))) * (x^(p - 1)) * exp(- (a * x + b / x) / 2)
    else
        return zero(T)
    end
end

function logpdf{T<:Real}(d::GeneralizedInverseGaussian{T}, x::Real)
    if x > 0
        a, b, p = params(d)
        return log(sqrt(a) / sqrt(b)) + Calculus.derivative(log(besselk(p, sqrt(a * b))) # How to get the derivative? Maybe log(pdf(x))
    else
        return -T(Inf)
    end
end

function cdf{T<:Real}(d::GeneralizedInverseGaussian{T}, x::Real)
  if x > 0
      a, b, p = params(d)
        return # (5) in Lemonte & Cordeiro 2011 Statistics & Probability Letters 81:506–517
    else
        return zero(T)
    end
end

@quantile_newton GeneralizedInverseGaussian

#### Sampling

# rand method from:
# Hörmann, W. & J. Leydold. (2014). Generating generalized inverse Gaussian random variates.
# J. Stat. Comput. 24: 547–557. doi:10.1007/s11222-013-9387-3

function rand(d::GeneralizedInverseGaussian)
    a, b, p = params(d)
    # Algorithm 1
    # Algorithm 2
    # Algorithm 3
end

Edit by ararslan: Wrap code in markdown code block

ararslan commented 7 years ago

Okay, thanks for your hard work here and sorry for the confusion!

nsgrantham commented 7 years ago

Thanks for kicking this off, @patwa67.

So far, I've defined the GeneralizedInverseGaussian with its initial methods (f14790479b7c8e5b50cddbdcd6bf6625f0067bc7) and implemented the sampling algorithms found in Hörmann & Leydold (2014) (c0bea481d4345bb426d9095d0a989be5aa9c7b81).

Next up is defining the cdf as given by eq. (5) in Lemonte & Cordeiro (2011).

What other methods ought to be included before I make a pull request?

JuliaStats / Distributions.jl

Generalized inverse Gaussian #554