Draw random numbers from a custom probability density function in Julia

I would like to draw random numbers from a modified exponential distribution:

p(x) = C * a * Exp[-(a*x)^b] with C=1/Gamma[1 + 1/b] for normalization.

How can I do this in julia? Unfortunately I have only little experience with Julia and no experiences with creating custom random numbers. I would be very grateful for any help.

Solution

If I'm not mistaken, that is a p-Generalized Gaussian distribution, which has a rather efficient implementation in Distributions.jl:

using Distributions
mu = 0       # your location parameter
alpha = 1/a  # your scale parameter
beta = b     # your shape parameter
p = PGeneralizedGaussian(mu, alpha, beta)

Using the Distributions.jl API for univariate distributions, you can sample from this distribution by passing it to the exported rand method. Here is an example of how to sample five independent scalars from a PGeneralizedGaussian distribution with mu = 0, alpha = 1/2 and beta = 3:

julia> p = PGeneralizedGaussian(0, 1/2, 3);

julia> rand(p, 5)
5-element Vector{Float64}:
  0.2835117212764108
 -0.023160728370422268
  0.3075395764050027
 -0.19233721955795835
  0.21256694763885342

If you want to try to implement the distribution yourself, which I do not recommend unless you are doing this as an exercise in programming math in Julia, you need to define a type that holds the static parameters of the distribution (in your case, just the shape parameter and the scale parameter), then define and export the methods listed here to extend the Distributions.jl API to your custom distribution. In particular, you need to define:

struct MyDistribution <: ContinuousUnivariateDistribution
  # ... your distribution parameters here
end

rand(::AbstractRNG, d::MyDistribution) # sample a value from d
sampler(d::MyDistribution)             # return d or something that can sample from d more efficiently
logpdf(d::MyDistribution, x::Real)     # compute the log of the pdf of d at x
cdf(d::MyDistribution, x::Real)        # compute the cdf of d at x
quantile(d::MyDistribution, q::Real)   # compute the qth quantile of d
minimum(d::MyDistribution)             # return the minimum value supported by d
maximum(d::MyDistribution)             # return the maximum value supported by d
insupport(d::MyDistribution, x::Real)  # query whether x is supported by d

The documentation of the package is very good, so it's an excellent way to get your feet wet if you are trying to learn Julia.