diff --git a/README.md b/README.md
index af67ef3..e9c8dc5 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,28 @@
# GeneralizedBeta1Distribution
+
+
[](https://mattiasvillani.github.io/GeneralizedBeta1Distribution.jl/stable/)
[](https://mattiasvillani.github.io/GeneralizedBeta1Distribution.jl/dev/)
[](https://github.com/mattiasvillani/GeneralizedBeta1Distribution.jl/actions/workflows/CI.yml?query=branch%3Amain)
[](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).
If `Y ~ Beta(α,β)` then `X = δ*Y^(1/γ)` follows `X ~ GeneralizedBeta1(α,β,γ,δ)`.
The implementation follows the Distributions.jl interface.
+
+## Installation
+Install from the Julia package manager (via Github) by typing `]` in the Julia REPL:
+```
+] add git@github.com: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
+```
\ No newline at end of file
diff --git a/src/GeneralizedBeta1.jl b/src/GeneralizedBeta1.jl
new file mode 100644
index 0000000..5a1f79c
--- /dev/null
+++ b/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))
+
diff --git a/src/GeneralizedBeta1Distribution.jl b/src/GeneralizedBeta1Distribution.jl
index 55bbd40..107d49d 100644
--- a/src/GeneralizedBeta1Distribution.jl
+++ b/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
diff --git a/test/GeneralizedBeta1Distribution.jl b/test/GeneralizedBeta1Distribution.jl
new file mode 100644
index 0000000..3df6ddd
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 450781b..db80a0d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,3 @@
-using GeneralizedBeta1Distribution
using Test
-@testset "GeneralizedBeta1Distribution.jl" begin
- # Write your tests here.
-end
+include("GeneralizedBeta1Distribution.jl")
\ No newline at end of file