first commit

README.md

@@ -1,6 +1,28 @@
 # GeneralizedBeta1Distribution
 
+
+
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://mattiasvillani.github.io/GeneralizedBeta1Distribution.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://mattiasvillani.github.io/GeneralizedBeta1Distribution.jl/dev/)
 [![Build Status](https://github.com/mattiasvillani/GeneralizedBeta1Distribution.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/mattiasvillani/GeneralizedBeta1Distribution.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/mattiasvillani/GeneralizedBeta1Distribution.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/mattiasvillani/GeneralizedBeta1Distribution.jl)
+
+## Description
+
+Julia implementation of the [Generalized Beta1 distribution (GB1)](https://en.wikipedia.org/wiki/Generalized_logistic_distribution#Type_IV).<br>If `Y ~ Beta(α,β)` then `X = δ*Y^(1/γ)` follows `X ~ GeneralizedBeta1(α,β,γ,δ)`.<br>The implementation follows the Distributions.jl interface.
+
+## Installation
+Install from the Julia package manager (via Github) by typing `]` in the Julia REPL:
+```
+] add [email protected]:mattiasvillani/GeneralizedBeta1Distribution.jl.git
+```
+
+## Example
+```julia
+using GeneralizedBeta1Distribution
+d = GeneralizedBeta1(1/2, 1/2, 2, 1)
+mean(d)    
+rand(d, 10)                                 # 10 random draws
+pdf(d, 0.5)                                 # pdf at 0.5
+cdf(d, 0.5)                                 # cdf at 0.5
+```

src/GeneralizedBeta1.jl

@@ -0,0 +1,142 @@
+import Distributions: logpdf, pdf, cdf, quantile, std, mean, skewness, @check_args, @distr_support
+import Base.rand
+import Base.length
+
+""" 
+    GeneralizedBeta1(α, β, γ, δ) 
+
+The Generalized Beta distribution of the first kind GB1 is a four-parameter extension of standard Beta(α, β) distribution.
+
+```math
+f(x) = \\frac{|γ|}{\delta^{\\alpha\\gamma}\\mathrm{B}(\\alpha,\\beta)} 
+       x^{\\alpha\\gamma - 1}(1 - (x/ \\delta)^ \\gamma)^{\\beta - 1}, \\quad 0 < x < 1
+```
+ 
+The distribution is defined by:
+
+If Y ∼ Beta(α, β), then X = (Y/δ)^γ ∼ GeneralizedBeta1(α, β, γ, δ).
+
+
+```julia
+d = GeneralizedBeta1(2, 3, 1.5, 1)   
+
+params(d)          # Get the parameters, i.e. α, β, γ and δ
+pdf(d, 0.5)        # Probability density function at x = 0.5
+mean(d)            # Mean 
+```
+
+External links
+
+* [Generalized Beta distribution of the first kind GB1 on Wikipedia](https://en.wikipedia.org/wiki/Generalized_beta_distribution#Generalized_beta_of_first_kind_(GB1))
+
+# Examples
+```julia-repl
+julia> d = GeneralizedBeta1(2, 3, 1.5, 1)
+julia> rand(d, 4)'
+1×4 adjoint(::Vector{Float64}) with eltype Float64:
+ 1.00851  0.640297  0.566234  2.16941
+julia> pdf(d, 1)
+julia> cdf(d, 1)
+```
+""" 
+struct GeneralizedBeta1{T<:Real} <: ContinuousUnivariateDistribution
+    α::T
+    β::T
+    γ::T
+    δ::T
+    GeneralizedBeta1{T}(α::T, β::T, γ::T, δ::T) where {T} = new{T}(α, β, γ, δ)
+end
+
+function GeneralizedBeta1(α::T, β::T, γ::T, δ::T; check_args::Bool=true) where {T <: Real}
+    @check_args GeneralizedBeta1 (α, α > zero(α)) 
+    @check_args GeneralizedBeta1 (β, β > zero(β))
+    @check_args GeneralizedBeta1 (δ, β > zero(δ))
+    return GeneralizedBeta1{T}(α, β, γ, δ)
+end
+
+GeneralizedBeta1(α::Real, β::Real, γ::Real, δ::Real; check_args::Bool=true) = 
+    GeneralizedBeta1(promote(α, β, γ, δ)...; check_args = check_args)
+GeneralizedBeta1(α::Integer, β::Integer, γ::Integer, δ::Integer; check_args::Bool=true) 
+    = GeneralizedBeta1(float(α), float(β), float(γ), float(δ); check_args = check_args)
+
+@distr_support GeneralizedBeta1 0 δ^γ
+
+Base.eltype(::Type{GeneralizedBeta1{T}}) where {T} = T
+
+#### Conversions
+convert(::Type{GeneralizedBeta1{T}}, α::S, β::S, γ::S, δ::S) where {T <: Real, S <: Real} = 
+GeneralizedBeta1(T(α), T(β), T(γ), T(δ))
+Base.convert(::Type{GeneralizedBeta1{T}}, d::GeneralizedBeta1) where {T<:Real} = 
+GeneralizedBeta1{T}(T(d.α), T(d.β), T(d.γ), T(d.δ))
+Base.convert(::Type{GeneralizedBeta1{T}}, d::GeneralizedBeta1{T}) where {T<:Real} = d
+
+#### Parameters
+params(d::GeneralizedBeta1) = (d.α, d.β, d.γ, d.δ)
+partype(::GeneralizedBeta1{T}) where {T} = T
+
+""" 
+    rand(d::GeneralizedBeta1[, n::Integer])
+
+Draw `n` random numbers from the GeneralizedBeta1 distribution `d`. 
+""" 
+function Base.rand(rng::Random.AbstractRNG, d::GeneralizedBeta1)
+    y = rand(rng, Beta(d.α, d.β))
+    return d.δ*y^(1/d.γ)
+end
+
+""" 
+    pdf(d::GeneralizedBeta1, x::Real) 
+
+Compute the pdf of the GeneralizedBeta1 distribution `d` at `x`. 
+""" 
+pdf(d::GeneralizedBeta1, x::Real) = ( abs(d.γ)/(d.δ^(d.α*d.γ)*beta(d.α,d.β)) )*
+    x^(d.α*d.γ - 1)*(1 - (x/d.δ)^d.γ)^(d.β - 1)
+
+
+""" 
+    logpdf(d::GeneralizedBeta1, x::Real) 
+
+Compute the logpdf of the GeneralizedBeta1 distribution `d` at `x`. 
+""" 
+logpdf(d::GeneralizedBeta1, x::Real) = log(abs(d.γ))-(d.α*d.γ)*log(d.δ) - logbeta(d.α, d.β) 
+    + (d.α*d.γ - 1)*log(x) + (d.β - 1)*log(1 - (x/d.δ)^d.γ)
+
+
+""" 
+    cdf(d::GeneralizedBeta1, x::Real) 
+
+Compute the cdf of the GeneralizedBeta1 `d` at `x`. 
+""" 
+cdf(d::GeneralizedBeta1, x::Real) = beta_inc(d.α, d.β, (x/d.δ)^d.γ)[1]
+
+
+""" 
+quantile(d::GeneralizedBeta1, p::Real) 
+
+Compute the `p`-quantile of the GeneralizedBeta1 distribution `d`. 
+""" 
+quantile(d::GeneralizedBeta1, p::Real) = d.δ * quantile(Beta(d.α, d.β), p)^(1/d.γ)
+
+""" 
+mean(d::GeneralizedBeta1) 
+
+Compute the mean of the GeneralizedBeta1 distribution `d`. 
+""" 
+mean(d::GeneralizedBeta1) = d.δ * (Beta(d.α + 1/d.γ, d.β)/Beta(d.α, d.β))
+
+
+""" 
+var(d::GeneralizedBeta1) 
+
+Compute the variance of the GeneralizedBeta1 distribution `d`. 
+""" 
+var(d::GeneralizedBeta1) = (d.δ^2) * (Beta(d.α + 2/d.γ, d.β)/Beta(d.α, d.β)) - mean(d)^2
+
+
+""" 
+    std(d::GeneralizedBeta1) 
+    
+Compute the standard deviation of the GeneralizedBeta1 distribution `d`. 
+""" 
+std(d::GeneralizedBeta1) = sqrt(var(d))
+

src/GeneralizedBeta1Distribution.jl

@@ -1,5 +1,17 @@
 module GeneralizedBeta1Distribution
 
-# Write your package code here.
+import Base: rand
+using Random
+import Random: default_rng, rand!
+import Statistics: mean, median, quantile, std, var
+import StatsBase: mode, params, params!
+using Distributions, SpecialFunctions
 
-end
+
+include("GeneralizedBeta1.jl")
+
+export mean, median, quantile, std, var, mode, params, params!, rand, rand!, skewness
+export pdf, cdf, logpdf, logcdf
+export GeneralizedBeta1
+
+end # End module

test/GeneralizedBeta1Distribution.jl

@@ -0,0 +1,33 @@
+using GeneralizedBeta1Distribution
+using Distributions: mean, std, var, pdf, logpdf, cdf, quantile
+
+@testset "GeneralizedBeta1Tests.jl" begin
+
+    d = GeneralizedBeta1(1/2,1/2,1/2,2)
+    @test cdf.(d, quantile.(d, 0.1:0.1:0.9)) ≈ 0.1:0.1:0.9
+
+    d = GeneralizedBeta1(3/2,3/2,2,2)
+    @test cdf.(d, quantile.(d, 0.1:0.1:0.9)) ≈ 0.1:0.1:0.9
+
+    @test pdf(d, 1) ≈ exp(logpdf(d, 1))
+
+    α = rand()
+    β = rand()
+    ds1 = GeneralizedBeta1(α,β,1,1) # should be the standard Beta(α,β)
+    ds2 = Beta(α,β) 
+    @test pdf(ds1, 0.5) ≈ pdf(ds2, 0.5) 
+    @test cdf(ds1, 0.5) ≈ cdf(ds2, 0.5) 
+
+    @test var(d) ≈ std(d).^2
+
+    @test length(rand(d, 4)) == 4
+
+    r = rand()
+    params(GeneralizedBeta1(2*r, r, 2, 2)) == (2*r, r, 2, 2)
+    @test mean(ds1) ≈ mean(ds2) ≈ α/(α+β)
+    @test var(ds1) ≈ var(ds2) ≈ (α*β)/((α+β)^2*(α+β+1))
+
+    @test GeneralizedBeta1(1, 2, 2, 2) == GeneralizedBeta1(1.0, 2.0, 2.0, 2.0)
+    @test GeneralizedBeta1(1, 2.0, 2, 2) == GeneralizedBeta1(1.0, 2.0, 2.0, 2.0)
+
+end

test/runtests.jl

@@ -1,6 +1,3 @@
-using GeneralizedBeta1Distribution
 using Test
 
-@testset "GeneralizedBeta1Distribution.jl" begin
-    # Write your tests here.
-end
+include("GeneralizedBeta1Distribution.jl")