User-defined distribution function example

question
#1

I am new to Julia, and I want to define my own distribution function. I found this post and followed the documentation, but I did not understand where do I put in the actual function (namely f(x)) of the distribution. Can someone help me with an example of say defining a gaussian distribution?

#2

You probably just need to learn how to create functions in general, as explained here. And then use that function for whatever you need, without worrying about the Distributions package. The idea of creating a sampler and all that is a little more involved than just a regular function. Do you know the equation you want to implement?

#3

@alejandromerchan my distribution is x\exp(-x) for x\in (0, \infty)
If I look at the example code from the post

type MyDist <: ContinuousUnivariateDistribution
    mu::Float64
    sigma::Float64
    MyDist(mu, sigma) = new(Float64(mu), Float64(sigma))
end

pdf(d::MyDist, x::Float64) = pdf(Normal(d.mu, 2*d.sigma), x)

I do not understand where should I introduce my function.
In this code the Normal distribution is used, but how to insert a custom function like f(x)=x*exp(-x) ?

#4

I’m not sure I understand the question, but maybe you just want

pdf(d::MyDist, x::Float64) = x * exp(-x)
#5

@jbrea I cannot even make the first part work, where the MyDist definition is.
Could you give an example of how to define a distribution with the function f(x)=xe^{-x}?

#6

Here is a full example

julia> using Distributions

julia> struct MyDist <: ContinuousUnivariateDistribution end

julia> Distributions.pdf(::MyDist, x) = x * exp(-x)
pdf (generic function with 73 methods)

julia> pdf(MyDist(), 3)
0.14936120510359183

The example you posted is quite old julia code; type was renamed to stuct at some point.

2 Likes
#7

@jbrea thanks!
But how do I sample from it?
I tried

struct MyDist <: ContinuousUnivariateDistribution end
Distributions.pdf(::MyDist, x) = x * exp(-x)
d = MyDist()
rand(d)

But this doesn’t work.
I thought I need to create a sampler according to the documentation, but I couldn’t figure out how to do it.
Could you add this to the example?

#8

As specified in the documentation you referred to, you need also to define

function Distributions.rand(::MyDist)
...
end

I guess inverse transform sampling could be useful.

#9

You can also import the desired functions if you don’t want to qualify with Distributions. For example,

import Distributions: rand,pdf,cdf

struct mydist <: ContinuousUnivariateDistribution
    ...
end 

function rand(dist::mydist)
    ...
end

...
#10

@Christopher_Fisher, @jbrea, my problem is that I’m not sure what should go into the .... I need an example of some implementation of such a rand function.

#11

The package Distributions.jl contains examples for all implemented distributions. You may use, e.g. @edit rand(Exponential(1.)) to see the implemention for the Exponential distribution or just browse the source code of the package.

#12

Btw. the distribution with pdf(x) = x * exp(-x) is known as the gamma distribution with shape α = 2 and rate β = 1 and it is already implemented as Gamma(2, 1).

#13

@jbrea, in the code I see

#### Sampling
rand(d::Exponential) = rand(GLOBAL_RNG, d)
rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))

I did not understand how can I use this.
BTW, what does the d before the :: stands for?

I probably need something like

function rand(d::MyDist)
    uni = Uniform()
    x = rand(uni)
    # do something with x, to get my distribution. e.g. x^2
end

I tried implementing the inverse transform sampling but here the CDF is 1- e^{-x}(x+1) which cannot be inverted.

Using the Gamma function is a nice solution, but I eventually want to add parameters to the function, which do not conform to the Gamma function.

#14

The argument. You should probably start with giving the manual thorough reading.

You can always find the root iteratively, using Newton’s method (usually works, and you can prove convergence) or bisection (worst case).

Providing all the methods is an investment that may not be needed for custom, one-off distributions, but this depends on the context.

#15

You may also want to have a look at wikipedia or this link, brought up by a simple web search.