Statistics of Cumulative sum of Poisson Process

Question

Hi I have posted this on the stack overflow but I think it is more appropriate to post here.

I'm trying to understand a link here between the Poisson and Exponential distributions. Here I have drawn values from a Poisson process and computed the cumulative sum

import numpy as np
import matplotlib.pyplot as plt

#Intensity of process
lam = 16
bins = 100
length_vec = 10000
poisson = np.random.poisson(lam,length_vec)
z_poisson = np.cumsum(poisson)
print('mean difference poisson cumsum = ',np.mean(np.diff(z_poisson)))
print('mean difference between poisson vals =',np.mean(np.abs(np.diff(poisson))))

#Histogram Poisson Process Vector
plt.hist(poisson, bins=bins)
plt.title('Poisson Process')
plt.xlabel('z')
plt.ylabel('count')

plt.show()
plt.close()

#Histogram Poisson Process Cumulative Sum Difference
plt.hist(np.diff(z_poisson), bins=bins)
plt.title('Poisson Process Interval Difference')
plt.xlabel('z')
plt.ylabel('count')

plt.show()

The first plot is a histogram of the Poisson process

and the second plot is a histogram of the difference of the elements of the cumulative sum of the Poisson process

If I take the mean of the difference between the cumulative sum elements vector, that value is

print('mean difference poisson cumsum = ',np.mean(np.diff(z_poisson)))
mean difference poisson cumsum = 16.000300030003

which is approximately almost equal to the intensity of the Poisson process. Does this mean that...

A) The increments of my cumulative sum z_{j+1}-z_j are Poisson distributed

B) The intensity of the process governs the size of the increments? I thought that the increments would be of size 1/lam ?

Looking for some general advice or reference to help me with this problem...can delete if not the correct stack. Thanks in advance.

Answer 1

That is not a Poisson process.

What you seem to have done is generate a sequence using the Poisson distribution, taken the cumulative sum, and then taken the differences of that cumulative sum; that should take you straight back to the original sequence (possibly missing the first value) so with the same distribution.

You can simulate a Poisson process by taking the cumulative sum of exponentially distributed random variables. If you then count the number of values in an interval, that should a random variable with a Poisson distribution, with mean equal to the rate of the exponentially distributed random variables multiplied by the width of the interval.

The following R code simulates that, with your rate of $16$ and an interval width of $1$, and also shows in red the theoretical Poisson distribution: the visual match is close given this is a simulation.

set.seed(2020) # for replication
rate <- 16
cases <- 100000 * rate
intervalwidth <- 1

process <- cumsum(rexp(cases, rate))                      # Poisson process
interval <- ceiling(process / intervalwidth)              # round to give a time interval
interval <- interval[interval != max(interval)]           # remove final part interval
count <- table(factor(interval, levels=1:max(interval)))  # count in interval
distribution <- table(factor(count, levels=0:max(count))) # table of counts
plot(distribution)                                        # graph of table of counts
points(0:max(count), dpois(0:max(count), rate * intervalwidth) * 
                     max(interval), col="red")            # theoretical distribution