How to add/subtract probability distributions?

[CentraCep]

As part of a mini-project where I am numerically solving a linear differential equation, I have to subtract a probability distribution from another distribution. Is there a way to do this in Julia? When I try:

a = Chi(3) - Uniform(0,1)

There is no method set up for this:

MethodError: no method matching -(::Chi{Float64}, ::Uniform{Float64})
Closest candidates are:
  -(::UnivariateDistribution, ::Real) at C:\Users\Acer\.julia\packages\Distributions\Fl5RM\src\univariate\locationscale.jl:139
  -(::ChainRulesCore.AbstractThunk, ::Any) at C:\Users\Acer\.julia\packages\ChainRulesCore\sHMAp\src\tangent_types\thunks.jl:30
  -(::ChainRulesCore.ZeroTangent, ::Any) at C:\Users\Acer\.julia\packages\ChainRulesCore\sHMAp\src\tangent_arithmetic.jl:101
  ...

[nilshg]

What is the result of Chi(3) - Uniform(0, 1)?

[juliohm]

You can take a look at the compositional data analysis literature. A PDF is an example of compositional data, and there is a space called the Aitchison space with well-defined vector addition, scalar multiplication and inner product:

https://github.com/JuliaEarth/CoDa.jl

The subtraction u - v is defined as usual as the addition of u with -v. You can see if this definition of subtraction is the one you are after.

[nilshg]

What does the result mean?

I had implicitly assumed that OP is talking about a convolution of probability distributions, although that of course only works in closed form for a few specific examples.

I guess the distribution of z = x - y with x \sim \chi(3) and y \sim U(0, 1) is this:

julia> using Distributions, Plots

julia> z = rand(Chi(3), 1_000_000) .- rand(Uniform(0, 1), 1_000_000);

julia> plot(x -> pdf(Chi(3), x), label = "Chi(3)")

julia> plot!(x -> pdf(Uniform(0, 1), x), label = "Uniform(0, 1)")

julia> histogram!(z, normalize = true, lw = 0.2, alpha = 0.5, label = "Chi(3) - Uniform(0, 1)")

[Tomas_Pevny]

Would this repo
https://github.com/cscherrer/Measures.jl
be of help?

[rafael.guerra]

Linking related thread.

[jonniedie]

If your plan is to eventually draw some fixed number of samples from the resulting distribution, you might want to check out MonteCarloMeasurements.jl. Doing it this way allows you to deal with complicated things like dependent random variables and operations that don’t have closed-form solutions. Take, for example, two uniformly distributed variables a∼U(1,2) and b∼U(1,2). Even though the two are identically distributed, a+b != a+a because a and a are (obviously) not independent.

julia> using Plots, MonteCarloMeasurements

julia> a = 1..2; b = 1..2
Particles{Float64, 2000}
 1.5 ± 0.289

julia> default(alpha=0.5, nbins=20, normalize=true)

julia> plot(a+a, label="a+a"); plot!(a+b, label="a+b")

[rafael.guerra]

@jonniedie, your example is splendid.
Please consider adding it to the Seven Lines of Julia collection.

[lrnv]

If this is convolution OP is after, then the right way to do it is probably via moment generating functions.

Of course it depends on what OP wants to do with the resulting distribution : sampling is easy. Computing the mgf is easy. Getting the PDF and/or CDF might pass through the Gaver-Stehfest algorithm (very fast convergence). Getting the quantile function might be easier by simulation dependeing on the needed precision.

There are specialised and efficient algorithm taylored for exactly what you need to do with the distribution, you just need to tell us.

What does OP want ?

[CentraCep]

@jonniedie @lrnv @juliohm @nilshg Sure, I’ll happily share more details about the nature of my task and hopefully it will provide sufficient information for someone to prescribe a solution.

I want to numerically solve the following linear differential equation in plasma physics that calculates the distribution in a collisional plasma, known as the BGK equation:

s denotes the particle species (e.g. protons, electrons, etc.), v_ss denotes the collisional frequency (this will be a constant, let’s say 1 for now), and F denotes a Maxwellian distribution (which is mathematically a chi distribution) that remains constant through time. The idea is that first, you can discretize this equation like I have done here explicitly:

v_col = 1

dt = 10^(-5) # timestep
t = range(0, 1.0, step=dt)

f_s = Vector{Uniform{Float64}}(undef, length(t)+1)
f_s[1] = Uniform(0,1)
F_s = Chi(3)

for it in 1:length(t)
    f_s[it + 1] = dt*v_col*(F_s - f_s[it]) + f_s[it]  
end

Then, you pick a distribution for f_s(0) (I am trying a uniform distribution first) and as you timestep through this equation, eventually what should happen is that the distribution f_s should relax into a Maxwellian F_s, which is what I am trying to show.

Before choosing this however, it is clear that two distributions have to be added or subtracted (F_s - f_s[it]). At the moment, when running this snippet, this error pops up

MethodError: no method matching -(::Chi{Float64}, ::Uniform{Float64})
Closest candidates are:
  -(::UnivariateDistribution, ::Real) at C:\Users\Acer\.julia\packages\Distributions\Fl5RM\src\univariate\locationscale.jl:139
  -(::ChainRulesCore.AbstractThunk, ::Any) at C:\Users\Acer\.julia\packages\ChainRulesCore\sHMAp\src\tangent_types\thunks.jl:30
  -(::ChainRulesCore.ZeroTangent, ::Any) at C:\Users\Acer\.julia\packages\ChainRulesCore\sHMAp\src\tangent_arithmetic.jl:101
  ...

I was wondering if someone has the best solution for this case. Perhaps this includes extraneous physics information, but I hope this addresses @lrnv 's request.

[lrnv]

Either you are not precise enough either I do not know your field enough to help.

What is the sum of distributions ? I do not know how to sum two distributions (and Julia neither, as your errors shows).

Are F and f random variables, which distributions are supposed to be chi2 and uniform ? In which case you should instentiate them through maybe MonteCarloMeasurements.jl

Another way to go at it would Be to consider a convolutional-friendly approximation of the Chi and Uniform distributions, e.g. the Laguerre basis (which has closed form convolutional properties, which would give you a gradient iteration that is easy to work with as you can precompute the expression of the terms in the basis and only keep the convolutional part in the main loop).

[jonniedie]

I don’t think there are any distributions besides gaussian that will remain closed under repeated addition/subtraction like this. You can use MonteCarloMeasurements, but I don’t know enough about your problem to know if F_s needs to be drawn once at the beginning or at each time step. You can use MonteCarloMeasurements in both cases, but it will take much longer if you need to redraw F_s every time step.

using Distributions, MonteCarloMeasurements, Plots

v_col = 1

dt = 10^(-5) # timestep
t = range(0, 1.0, step=dt)

Here is the method with a single draw:

f_s0 = Particles(Uniform(0,1))
f_s = [f_s0]
F_s = Particles(Chi(3))

for it in 1:length(t)-1
    push!(f_s, dt*v_col*(F_s - f_s[it]) + f_s[it])
end

plot(t, f_s)

and the method with a new draw for each time step (but keep in mind, this is going to be really slow):

f_s0 = Particles(Uniform(0,1))
f_s = [f_s0]

for it in 1:length(t)-1
    push!(f_s, dt*v_col*(Particles(Chi(3)) - f_s[it]) + f_s[it])
end

plot!(t, f_s)

I’m not sure which of those looks like what you were expecting.

But all of this is much slower than you’d probably want it to be. You may want to consider using DifferentialEquations.jl for the adaptive time-stepping with an EnsembleProblem if you’re going the Monte Carlo route. You could also swap out Particles for StaticParticles for speed, but this will give you fewer traces.

[CentraCep]

For the plot, I would like the dependent variable to be the probability density and the independent variable to be the different outcomes (or in this case, different velocities).

So for example, initially, before timestepping, my distribution f_s and my Maxwellian F_s look like this:

# Initially comparing chosen distribution f_s to Maxwellian F_s

velᵤ = -0.1:0.005:4
f_s0 = pdf.(Uniform(0,1), velᵤ) # uniform distribution with area 1
plot(velᵤ, f_s0, label=L"f_{s}0", framestyle=:box)

velᵪ = 0:0.005:4
F_s = pdf.(Chi(3), velᵪ)  # chi distribution with area 1

plot!(velᵪ, F_s, linecolor = :orange, linestyle = :dash, label= L"F_s" * " (Maxwellian)")
xlabel!("velocity")
ylabel!("probability density")
xlims!(-0.1, 4)
ylims!(0, 1.5)

Then ideally, by the end of time, the blue line (f_s) should overlap the orange line (F_s).

[dlakelan]

@CentraCep

I want to numerically solve the following linear differential equation in plasma physics that calculates the distribution in a collisional plasma, known as the BGK equation:

s denotes the particle species (e.g. protons, electrons, etc.), v_ss denotes the collisional frequency (this will be a constant, let’s say 1 for now), and F denotes a Maxwellian distribution (which is mathematically a chi distribution) that remains constant through time. The idea is that first, you can discretize this equation like I have done here explicitly:

It’s interesting how ambiguous mathematical notation can be.

I’m going to go out on a limb here and say that this is supposed to be a description of point-wise behavior of the density function f(v). As such you could think of it as an infinite number of ordinary differential equations (one for each value of v)?

reading further in your comment I believe I’m correct. So the part you seem to be missing is that in addition to discretizing time, you also need to discretize “space” or in this case “velocity” and you have to evaluate the pdf of the distribution at each v value to form a bunch of individual differential equations one for each v range.

In addition to this problem, you have to ensure that the distributions maintain their “sum to 1” constraint. Because of the sum to 1 constraint it’s entirely plausible that the equation can’t actually hold pointwise for every point in the velocity domain, so you will have to decide what the equation even means, perhaps for example it holds on average in the sense of an integral across the velocity domain?

[rafael.guerra]

@CentraCep

differential equation in plasma physics

In case it helps, linking some Julia plasma physics resources: Julia Vlasov, Plasmas.jl and Plasma.jl.

[CentraCep]

Also, when I try running your snippets, I get:

BoundsError: attempt to access 100000-element StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64} at index [1:100001]

[jonniedie]

Oh, yeah f_s has to have the same length as the time vector for plotting. I thought I fixed that, but something may be off.

[CentraCep]

Figured out the issue, it was an off by 1 error with length(n).

Anyway, please refer to my reply from earlier with the plots of both distributions at time 0: I would like the dependent variable to be the probability density and the independent variable to be the different outcomes (or in this case, different velocities).

Also, F_s only needs to be drawn once.

[jonniedie]

Using my first method above (I went back and edited it to fix the off-by-one error), you get this:

Does that look close enough to what you were envisioning? It’s hard to see because it moves so quickly, but the initial distribution is the same uniform one you plotted above.

[CentraCep]

Yes, this was exactly what I was envisioning!

May I ask how you created this plot/animation in Julia?