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.

using Distributions
using Plots

The function below makes the simulation of a homogeneous Poisson process with given intensity lambda and time horizon tmax.

function HPoisson_T(tmax, lambda)
    t=0
    output = [] # intialize the vector of points to be returned
    while t <= tmax
        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)

Exercise 2 (GOF tests)

Part 1

Simulate several trajectories of a homogeneous Poisson process with rate $\lambda=1$ on $[0,T]$. Perform one of the 3 tests presented during the course (Uniform test, Exponential test, Donsker test) on these trajectories.

As we saw during the lecture, the Exponential test is the most powerful in order to detect clustering or repulsiveness.

Here is the code for the Exponential test. It involves standard Kolmogorov Smirnov test and so we use the package HypothesisTests.

using HypothesisTests

function ExponentialTest(points)
    isempty(points) && return 0
    S = diff([0;points])
    pvalue(ExactOneSampleKSTest(S, Exponential()))
end

ExponentialTest (generic function with 1 method)

Question 1

Compute the empirical type I error of these procedures while $T$ is varying from $10$ to $1000$.

Nsimu = 5e2
lambda = 1
T = 1e3
Tvec = LinRange(10,T,10)
reject = zeros(length(Tvec))
for _ in 1:Nsimu
    for k in eachindex(Tvec)
        points = HPoisson_T(Tvec[k], lambda)
        reject[k] += (ExponentialTest(points) < .05)/Nsimu
    end
end
scatter(Tvec, 100*reject, label = "", title = "Empirical type I error (λ="*string(lambda)*")", xlab="T", ylab="%", ylims=[0,10])
plot!([0,T], [5,5], label="5%")

Question 2

Compute the empirical type II error of these procedures with $\lambda\neq 1$. How does it vary with respect to $T$ ?

Nsimu = 2e2
lambda = 1.1
T = 1e3
Tvec = LinRange(10,T,10)
reject = zeros(length(Tvec))
for _ in 1:Nsimu
    for k in eachindex(Tvec)
        points = HPoisson_T(Tvec[k], lambda)
        reject[k] += (ExponentialTest(points) < .05)/Nsimu
    end
end
scatter(Tvec, 100*(1 .- reject), legend=false, title = "Empirical type II error (λ="*string(lambda)*")", xlab="T", ylab="%", ylims=[0,100])

The type II error is decreasing as a function of $T$. This is natural since the value of $T$ is related with the sample size $n$ in the standard statistics framework.

Part 2

Simulate several trajectories of the non homogeneous Poisson process described there with $T=\alpha=\beta=1$. After using the time change Theorem, perform one of the 3 tests presented during the course (Uniform test, Exponential test, Donsker test) on these trajectories.

The cumulative intensity of this model 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.

Here is the code to simulate such a process.

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 1 method)

Here is a custom version of the Beta incomplete function.

using SpecialFunctions

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

beta_inc_perso (generic function with 1 method)

Here is the code to perform the change of time.

function changetime(points, T, C, alpha, beta)
    return C*T .* beta_inc_perso.(alpha+1, beta+1, points ./ T) .* SpecialFunctions.beta(alpha+1, beta+1)
end

changetime (generic function with 1 method)

Question 1

Compute the empirical type I error of these procedures while $C$ is varying from $20$ to $1000$.

T=1; alpha=1; beta=1; alpha0=1; beta0=1
Nsimu = 1e3
C = 1e3
Cvec = LinRange(10,C,10)
reject = zeros(length(Cvec))
for _ in 1:Nsimu
    for k in eachindex(Cvec)
        points = NHPoisson_conditional(T, Cvec[k], alpha, beta)
        scaledpoints = changetime(points, T, Cvec[k], alpha0, beta0)
        reject[k] += (ExponentialTest(scaledpoints) < .05)/Nsimu
    end
end
scatter(Cvec, 100*reject, label="", title = "Empirical type I error (α="*string(alpha)*", β="*string(beta)*")", xlab="C", ylab="%", ylims=[0,10])
plot!([0,C], [5,5], label="5%")

Like in the first part, the test is exact.

Question 2

Compute the empirical type II error of these procedures when the values of $\alpha$ and $\beta$ are misspecified. How does it vary with respect to $C$ ?

T=1; alpha=1; beta=.8; alpha0=1; beta0=1
Nsimu = 1e3
C = 1e3
Cvec = LinRange(20,C,10)
reject = zeros(length(Cvec))
for _ in 1:Nsimu
    for k in eachindex(Cvec)
        points = NHPoisson_conditional(T, Cvec[k], alpha, beta)
        scaledpoints = changetime(points, T, Cvec[k], alpha0, beta0)
        reject[k] += (ExponentialTest(scaledpoints) < .05)/Nsimu
    end
end
scatter(Cvec, 100*(1 .- reject), label="", title = "Empirical type II error (α="*string(alpha)*", β="*string(beta)*")", xlab="C", ylab="%", ylims=[0,100])

T=1; alpha=1.3; beta=1; alpha0=1; beta0=1
Nsimu = 1e3
C = 1e3
Cvec = LinRange(20,C,10)
reject = zeros(length(Cvec))
for _ in 1:Nsimu
    for k in eachindex(Cvec)
        points = NHPoisson_conditional(T, Cvec[k], alpha, beta)
        scaledpoints = changetime(points, T, Cvec[k], alpha0, beta0)
        reject[k] += (ExponentialTest(scaledpoints) < .05)/Nsimu
    end
end
scatter(Cvec, 100*(1 .- reject), label="", title = "Empirical type II error (α="*string(alpha)*", β="*string(beta)*")", xlab="C", ylab="%", ylims=[0,100])

In comparison with the first part, the role of $T$ (merely the control of the sample size) is replaced by the role of $C$. With this comparison in mind, we observe here also that the type II error is a decreasing function of $C$.