# Poisson: Probability distribution of amount of fuel per day based on arrival rate of cars at gas station.

I am modeling a hydrogen refueling station with an average of 2.5 cars arriving per day. This can be modeled using a Poisson distribution.

Each car has the following (discrete) probability distribution for the amount of hydrogen that it refuels. This probability distribution is based on a histogram derived from historical data.

How can I create a probability distribution for the total amount of hydrogen tanked per day?

I have been thinking how I can "add" probability distributions but cannot find a sensible way to do this. I assume it has something to do with the convolution of a Poisson distribution with another distribution.

The purpose is to create a discrete probability distribution and use it as an input to the Stochastic Dual Dynamic Programming (SDDP) algorithm implemented in SDDP.jl in order to optimize an energy system. So, my desired outcome is something like this, but then of course with correct data (i.e. no uniform distribution).

• Welcome to MSE. Your question is phrased as an isolated problem, without any further information or context. This does not match many users' quality standards, so it may attract downvotes, or be closed. To prevent that, please edit the question. This will help you recognize and resolve the issues. Concretely: please provide context, and include your work and thoughts on the problem. These changes can help in formulating more appropriate answers. Commented Apr 3, 2023 at 9:54
• Various issues: (a) What have you tried? (b) Your normal distribution needs a standard deviation or variance if you want anything more than the mean. (c) Normal distributions have a positive probability of being negative. (d) the sum of a Poisson number of normal distributions is unlikely to have a simple distribution - what information do you want? (e) Simulation might be faster than analysis Commented Apr 3, 2023 at 9:59
• edited the question Commented Apr 3, 2023 at 10:56

The resulting distribution is a discrete compound Poisson distribution.

The python code below can create the distribution. Be aware that this code is not specific to the Poisson process but can take any distribution for the amount of cars per day.

import pandas as pd
import numpy as np

def poisson_plus_discrete(n_cars_dist, kg_dist, kg_step=1, verbose=False):
#lower bound of the amount of cars and kg must be zero!
#discrete compound poisson distribution
if sum(n_cars_dist) != 1:
print('ERROR: probabilities do not add up to 1.')
max_cars = len(n_cars_dist) - 1

cols = np.arange(0, max_cars*kg_step*(len(kg_dist)-1) + kg_step/2, kg_step).tolist()
condf = pd.DataFrame(columns=['prob']+cols)
condf.index.name = 'n_cars'
kg_row_0 = padding([1.0], len(condf.columns) - 1)
condf.loc[0] = np.insert(kg_row_0, 0, n_cars_dist[0])
kg_row_1 = padding(kg_dist, len(condf.columns) - 1)
condf.loc[1] = np.insert(kg_row_1, 0, n_cars_dist[1])
probabilities = n_cars_dist[0] * kg_row_0 + n_cars_dist[1] * kg_row_1

convolved = kg_dist
for i in range(2, max_cars + 1):
convolved = np.convolve(convolved, kg_dist) #new length is N + M - 1
kg_row = padding(convolved, len(condf.columns) - 1)
condf.loc[len(condf)] = np.insert(kg_row, 0, n_cars_dist[i])
probabilities = probabilities + n_cars_dist[i] * kg_row

if verbose:
print(condf)
print('\nprobabilities:', probabilities.tolist())
print('kg of hydrogen:', cols)

return probabilities

n_cars_dist = [0.3, 0.3, 0.4]
#n_cars_dist = [0.3, 0.3, 0.2, 0.1, 0.05, 0.05] # for the values [0,1,2]
kg_dist = [0.1, 0.6, 0.3] # for the values [0,1,2]
kg_step = 1
poisson_plus_discrete(n_cars_dist, kg_dist, kg_step, verbose=True)


What is described is a compound Poisson distribution.

It's straightforward to numerically generate a discrete probability distribution, by successive convolutions of the refueling distribution with itself and truncating the number of cars arriving per day at a relatively large number (sorry, it's in R, not Julia):

library(data.table)
set.seed(34836697)

# dummy refuel mass probabilities for 0.2:20 kg
m <- runif(50)
m <- m/sum(m)

lambda <- 2.5 # average cars/day

system.time({
g200_0 <- 1:50 # tanked amount/car in units of 200 g
# Truncate the number of cars/day. The probability of seeing more cars than
# this is less than machine precision.
maxcars <- qpois(.Machine$double.eps, lambda, 0) # initialize a list to hold data.tables of total daily tanked probabilities # for each of 0:maxcars number of cars in a day ldt <- vector("list", maxcars + 1L) # if no cars come the tanked amount is 0 with probability 1 ldt[[1]] <- data.table(g200 = 0L, p = 1, cars1 = 1L) # get the distribution of tanked amount conditioned on the total number of # cars arriving in a day for (i in 2:(maxcars + 1L)) { ldt[[i]] <- ldt[[i - 1L]][ ,.(p = c(outer(p, m)), g200 = c(outer(g200, g200_0, "+"))) ][, .(p = sum(p), cars1 = i), g200] } d <- dpois(0:maxcars, lambda) # total daily car visit probabilities # get the marginal distribution of mass tanked dt <- setnames( setorder(rbindlist(ldt)[, .(p = sum(p*d[cars1])), g200], g200)[ , g200 := g200/5 # convert to kg ], "g200", "kg" ) }) #> user system elapsed #> 0.07 0.00 0.08  dt contains the probabilities of daily tanked amounts in kg sum(dt$p)
#> [1] 1

dt[]
#>          kg            p
#>    1:   0.0 8.208500e-02
#>    2:   0.2 3.700003e-03
#>    3:   0.4 8.753526e-04
#>    4:   0.6 8.869384e-03
#>    5:   0.8 3.139474e-03
#>   ---
#> 1197: 239.2 2.296300e-48
#> 1198: 239.4 5.492619e-49
#> 1199: 239.6 1.250848e-49
#> 1200: 239.8 1.860561e-50
#> 1201: 240.0 3.030906e-51

• You are pointing me in the right direction. However, it turns out that it should be a discrete compound Poisson distribution en.wikipedia.org/wiki/Compound_Poisson_distribution. In the mixed Poisson distribution, the rate parameter is random, which is not the case in my question. I can't read R well but will post my Python solution below. Commented Apr 4, 2023 at 12:07
• You are correct. It looks like your python code is doing something very similar with successive convolutions. Commented Apr 4, 2023 at 12:34