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 Likes
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?
1 Like
@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) ?
I’m not sure I understand the question, but maybe you just want
pdf(d::MyDist, x::Float64) = x * exp(-x)
1 Like
@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}?
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.
3 Likes
@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?
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.
1 Like
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
...
@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.
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.
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).
@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.
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.