# How to add/subtract probability distributions?

**URL:** <https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673>\
**Category:** New to Julia\
**Tags:** physics, distributions\
**Created:** [February 2, 2022, 4:25pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673 "2022-02-02T16:25:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 2, 2022, 4:25pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/1 "2022-02-02T16:25:20Z")

</div>

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:

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

```

There is no method set up for this:

```julia
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
  ...

```

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [February 2, 2022, 5:10pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/2 "2022-02-02T17:10:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [February 2, 2022, 6:42pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/4 "2022-02-02T18:42:11Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [February 2, 2022, 6:47pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/5 "2022-02-02T18:47:51Z")

</div>

What does the result mean?

I had implicitly assumed that OP is talking about a [convolution of probability distributions](https://en.wikipedia.org/wiki/List_of_convolutions_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
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)")

```

![image](https://global.discourse-cdn.com/julialang/original/3X/7/9/79d276cfef5f00e2594be90c1318580490b9a8bd.png)

---

<div class="post-metadata">

**Author:** ![Tomas\_Pevny](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomas_pevny/32/25466_2.png) [@Tomas\_Pevny](https://discourse.julialang.org/u/Tomas_Pevny)\
**Post date:** [February 2, 2022, 6:49pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/6 "2022-02-02T18:49:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 2, 2022, 7:01pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/7 "2022-02-02T19:01:54Z")

</div>

Linking [related thread](https://discourse.julialang.org/t/distribution-of-sum-of-random-variables/69023).

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [February 2, 2022, 7:28pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/8 "2022-02-02T19:28:09Z")

</div>

If your plan is to eventually draw some fixed number of samples from the resulting distribution, you might want to check out [MonteCarloMeasurements.jl](https://github.com/baggepinnen/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
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")

```

![image](https://global.discourse-cdn.com/julialang/original/3X/9/d/9d582ed8cd8c12ee822d2ee734f9713925d1c3c3.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 2, 2022, 7:52pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/9 "2022-02-02T19:52:56Z")

</div>

@jonniedie, your example is splendid.  
Please consider adding it to the [_Seven Lines of Julia_](https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416) collection.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [February 3, 2022, 12:46am UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/10 "2022-02-03T00:46:49Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 4, 2022, 8:17pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/11 "2022-02-04T20:17:55Z")

</div>

@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:

![BGK](https://global.discourse-cdn.com/julialang/original/3X/c/1/c16b2ec229ac04d2589718360922da967e3ace95.png)

`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](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution)) that remains constant through time. The idea is that first, you can discretize this equation like I have done here explicitly:

```julia
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

```julia
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.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [February 5, 2022, 2:20am UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/12 "2022-02-05T02:20:20Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [February 5, 2022, 3:00am UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/13 "2022-02-05T03:00:52Z")

</div>

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.

```julia
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:

```julia
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):

```julia
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/7/47be1c28ade5991e5b542894a480d90905167b80.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 5, 2022, 6:41pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/14 "2022-02-05T18:41:40Z")

</div>

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:

```julia
# 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)

```

 ![Screenshot 2022-02-05 124046](https://global.discourse-cdn.com/julialang/original/3X/9/5/953477f05bfd6477c8dc7a3a7ef30f28bc887d6c.png)

Then ideally, by the end of time, the blue line (f\_s) should overlap the orange line (F\_s).

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [February 5, 2022, 7:40pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/15 "2022-02-05T19:40:08Z")

</div>

> [@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:
> 
> ![BGK](https://global.discourse-cdn.com/julialang/original/3X/c/1/c16b2ec229ac04d2589718360922da967e3ace95.png)
> 
> `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](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_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?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 6, 2022, 12:19am UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/16 "2022-02-06T00:19:26Z")

</div>

> [@CentraCep](#):
>
> differential equation in plasma physics

In case it helps, linking some Julia plasma physics resources: [Julia Vlasov](https://juliavlasov.github.io/), [Plasmas.jl](https://github.com/djsegal/Plasmas.jl) and [Plasma.jl](https://github.com/killah-t-cell/Plasma.jl).

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 6, 2022, 8:13pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/17 "2022-02-06T20:13:20Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [February 6, 2022, 8:28pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/18 "2022-02-06T20:28:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 6, 2022, 8:37pm UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/19 "2022-02-06T20:37:06Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jonniedie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jonniedie/32/12842_2.png) [@jonniedie](https://discourse.julialang.org/u/jonniedie)\
**Post date:** [February 7, 2022, 4:50am UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/20 "2022-02-07T04:50:59Z")

</div>

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

![probs](https://global.discourse-cdn.com/julialang/original/3X/3/a/3ab5d2f5e678dfa4e72f0b6bee35dbe10bcc8ce8.gif)

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.

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 7, 2022, 5:15am UTC](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673/21 "2022-02-07T05:15:11Z")

</div>

Yes, this was exactly what I was envisioning! 🙂

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

[Next page](https://discourse.julialang.org/t/how-to-add-subtract-probability-distributions/75673.md?page=2)