Preamble¶

The code is written in Julia programming language. However, I think that the code is readable without any knowledge in Julia and can be easily adapted to Python or R for instance.

Here is the preamble to load the Distributions and Plots packages.

In [ ]:
using Distributions
using Plots
using SpecialFunctions

Here is a function used in the following in order to get the x and y coordinates to plot the graph of a counting process.

In [ ]:
function xy_countingprocess(tmin, tmax, points)
    n = length(points)
    x = points
    y = collect(1:length(points))
    pushfirst!(x, tmin); pushfirst!(y, 0) # adds the t=tmin values to the x and y coordinates
    push!(x, tmax); push!(y, n) # adds the t=tmax values to the x and y coordinates
    return (x,y)
end

xy_countingprocess (generic function with 1 method)

Exercise 5¶

Here are two algorithms based on the renewal structure of the homogeneous Poisson process. The first one uses a prescribed number $n$ of points to be simulated. The second one uses a prescribed time horizon $T$ for the simulation.

In [ ]:
function HPoisson_n(n, lambda)
    if n==0  # used to ensure that the output is a Vector (even if n=0 or 1)
        return []
    end
    if n==1
        return [rand(Exponential(1/lambda))]
    end
    S = rand(Exponential(1/lambda), n)
    output = cumsum(S)
    return output
end

function HPoisson_T(T, lambda)
    t=0
    output = [] # intialize the vector of points to be returned
    while t <= T
        S = rand(Exponential(1/lambda))
        t += S # updates the time value
        push!(output, t) # appends the value of t to the output vector
    end
    pop!(output) # removes the last element (which is strictly large than T)
    return output
end

HPoisson_T (generic function with 1 method)

Question 1¶

Here is the code to generate a graph of the counting processes (CP) between $t=0$ and $t=30$. For each process, its cumulative intensity (CI) is given in transparency.

In [ ]:
T=30.

colors = distinguishable_colors(3, RGB(0,0,1)) # Used to get 3 distinguishable colors

lambda = .5
points = HPoisson_T(T, lambda)
(x,y) = xy_countingprocess(0., T, points)
plot(x, y, linetype = :steppost, color = colors[1], label = "CP, λ=.5")
plot!([0,T], [0,lambda*T], color = colors[1], alpha = .3, label = "CI, λ=.5")

lambda = 1.
points = HPoisson_T(T, lambda)
(x,y) = xy_countingprocess(0., T, points)
plot!(x, y, linetype = :steppost, color = colors[2], label = "CP, λ=1.")
plot!([0,T], [0,lambda*T], color = colors[2], alpha = .3, label = "CI, λ=1.")

lambda = 2.
points = HPoisson_T(T, lambda)
(x,y) = xy_countingprocess(0., T, points)
plot!(x, y, linetype = :steppost, color = colors[3], label = "CP, λ=2.")
plot!([0,T], [0,lambda*T], color = colors[3], alpha = .3, label = "CI, λ=2.")

One can see the expected behavior that, the larger the intensity is, the larger the number of points is. Moreover, one can see that the counting process is close to the graph of the cumulative intensity (the 3 straight lines here).

Question 2¶

In this question, we estimate some expectation values via Monte Carlo method. Let us fix $\lambda$ and the number of simulations used for Monte Carlo.

In [ ]:
lambda = 2.
Nsimu = 1e3

1000.0

a) Estimation of $E[T_2]$.¶

In [ ]:
MCmean = 0
for k in 1:Nsimu
    points = HPoisson_n(2,lambda)
    T2 = points[2]
    MCmean += T2/Nsimu
end
MCmean

0.9750693993897449

b) Estimation of $E[N_{10}]$.¶

In [ ]:
MCmean = 0
for k in 1:Nsimu
    points = HPoisson_T(10,lambda)
    N10 = length(points)
    MCmean += N10/Nsimu
end
MCmean

20.071999999999907

Exercise 6¶

Here is an algorithm based on Lewis-Shedler thinning method.

In [ ]:
function NHPoisson_Thinning(T, lambda_function, lambda_bound)
    points_bound = HPoisson_T(T, lambda_bound) # These are the points that (a priori) belong to the output
    output = []
    n = length(points_bound)
    U = rand(Uniform(0, lambda_bound), n) # Draw a uniform rv for each (a priori) point
    for i in 1:n
        if U[i] < lambda_function(points_bound[i])
            push!(output, points_bound[i]) # if the uniform is less than the intensity we keep the point
        end
    end
    return output
end

NHPoisson_Thinning (generic function with 1 method)
In [ ]:
T = 30
lambda_function(t) = 2 .+ 2 .* sin(t)
lambda_bound = 4

cumulative_intensity(t) = 2 .* t .+ 2 .* (1 .- cos.(t))

points = NHPoisson_Thinning(T, lambda_function, lambda_bound)
(x,y) = xy_countingprocess(0., T, points)
plot(x, y, linetype = :steppost, color = :blue, label = "NH-PP")
xvec = LinRange(0, T, 1000)
plot!(xvec, cumulative_intensity(xvec), color = :blue, alpha = .3, label = "CI")

Exercise 7¶

Here is an algorithm based on the conditional structure of the Poisson process. For ease of presentation, it assumes that the cumulative intensity function $\Lambda$ and a simulation procedure for the $n$-sample with density $f$ are given.

In [ ]:
function NHPoisson_conditional(T, cumulative_intensity, simulation_procedure)
    n = rand(Poisson(cumulative_intensity(T)))
    output = simulation_procedure(n)
    sort!(output)
    return output
end

NHPoisson_conditional (generic function with 2 methods)

However, this general presentation is not usable in practice. Here is below an adaptation to the framework of Exercice 7.

We use the fact that $\Lambda(T) = \int_0^T C(t/T)^\alpha (1-t/T)^\beta dt = CT \int_0^1 x^\alpha (1-x)^\beta dx = CT \cdot \operatorname{B}(\alpha+1, \beta+1)$, where $\operatorname{B}$ is the bêta function.

Moreover, the density $f$ simplifies to $f(t) = \frac{(t/T)^\alpha (1-t/T)^\beta}{T \operatorname{B}(\alpha+1, \beta+1)} \mathbf{1}_{[0,T]}(t)$ which is the density of the random variable $T\cdot \xi$ where $\xi$ is a random variable with bêta distribution with parameters $\alpha+1$ and $\beta+1$.

In [ ]:
function NHPoisson_conditional(T, C, alpha, beta)
    Lambda = C*T*SpecialFunctions.beta(alpha+1, beta+1)
    n = rand(Poisson(Lambda))
    if n==0 # used to ensure that the output is a Vector (even if n=0 or 1)
        return []
    end
    if n==1
        return [T * rand(Beta(alpha+1, beta+1))]
    end
    output = T .* rand(Beta(alpha+1, beta+1), n)
    sort!(output) 
    return output
end

NHPoisson_conditional (generic function with 2 methods)
In [ ]:
Cvec = [1, 10]
T = 30.
alphavec = [.5, 1, 2]
betavec = [.5, 1, 2]

colors = distinguishable_colors(9)
i = 0
plot()
for (alpha, beta, C) in Iterators.product(alphavec, betavec, Cvec)
    points = NHPoisson_conditional(T, C, alpha, beta)
    (x,y) = xy_countingprocess(0., T, points)
    plot!(x, y, linetype = :steppost, color = colors[mod(i,9)+1], label = "C = "*string(C)*", α = "*string(alpha)*", β = "*string(beta))
    i += 1
end
display(plot!())

Of course, the $C$ parameter only controls the mean number of points, whereas $\alpha$ and $\beta$ control the shape of the intensity. In the plot above, there are two curves for each color: they are associated with the same parameters $(\alpha, \beta)$ and, the lower curve corresponds to $C=1$ and the upper curve corresponds to $C=10$.

Here are two plots to understand the behavior of the process with respect to $\alpha$ and $\beta$ (for $C=1$ for instance).

  1. the cumulative intensity,
  2. the intensity.

As in the previous, we could compare the previous plot (counting processes) with the the first plot below.

Remarks:

  • The cumulative intensity is $\Lambda(s) = \int_0^s C(t/T)^\alpha (1-t/T)^\beta dt = CT \int_0^{s/T} x^\alpha (1-x)^\beta dx = CT \cdot I_{s/T}(\alpha+1,\beta+1)\operatorname{B}(\alpha+1, \beta+1)$, where $I_x$ is the regularized incomplete bêta function.
  • $I_x$ (as a function of $x$) is a c.d.f. Hence, as appear in the expression above, the intensity is not normalized from one choice of $(\alpha,\beta)$ to the other. This was apparent in the previous plot and is very striking in the next one (the terminal value of the cumulative intensity depends on $(\alpha,\beta)$).

1. Cumulative intensity¶

In [ ]:
T = 30.
alphavec = [.5, 1, 2]
betavec = [.5, 1, 2]

function beta_inc_perso(alpha, beta, x) # the beta_inc function of the SpectialFunctions package returns a tuple which is annoying for our utilization
    return SpecialFunctions.beta_inc(alpha, beta, x)[1]
end

x = LinRange(0,T,100)
colors = distinguishable_colors(9)
i = 0
plot()
for (alpha, beta) in Iterators.product(alphavec, betavec)
    y = T .* beta_inc_perso.(alpha+1, beta+1, x/T) .* SpecialFunctions.beta(alpha+1, beta+1)
    plot!(x, y, color = colors[i+1], label ="α = "*string(alpha)*", β = "*string(beta))
    i += 1
end
display(plot!())

2. Intensity¶

In [ ]:
T = 30.
alphavec = [.5, 1, 2]
betavec = [.5, 1, 2]

x = LinRange(0,T,100)
colors = distinguishable_colors(9)
i = 0
plot()
for (alpha, beta) in Iterators.product(alphavec, betavec)
    y = (x/T).^alpha .* (1 .- x/T).^beta
    plot!(x, y, color = colors[i+1], label ="α = "*string(alpha)*", β = "*string(beta))
    i += 1
end
display(plot!())
In [ ]: