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

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 pointprocess_withintensity(tmin, tmax, points, intensity)
    npoints = length(points)
    scatter(points, zeros(npoints))
    tvec = LinRange(tmin, tmax, 1000)
    intensityvec = Float64[]
    for i in eachindex(tvec)
        append!(intensityvec, intensity(tvec[i], points[points .< tvec[i]]))
    end
    plot!(tvec, intensityvec, legend = false)
end

pointprocess_withintensity (generic function with 1 method)

Exercise 4 (Thinning simulation)¶

Use the idea of the thinning representation to write an algorithm which gives a simulation of a point process.

The function below gives the simulation of a point process via the thinning algorithm. It assumes that a bound for the intensity (at least up to some future time) is known.

The arguments intensity and intensity bound are functions of t and of the previous points.

  • intensity returns a scalar value,
  • intensitybound returns a scalar value and a future time after which the bound is not valid anymore.
In [ ]:
function thinning(tmax, intensity, intensitybound)
    t = 0
    points = Real[]
    while t < tmax
        (bound, tbound) = intensitybound(t, points)
        IEI = rand(Exponential(1/bound))
        t += IEI
        t > tmax && break
        if t > tbound
            t = tbound
            continue
        end
        truerate = intensity(t, points)
        truerate > bound && error("intensitybound must be an upperbound of the function intensity")
        truerate < 0 && error("intensity must be non negative")
        U = rand(Uniform(0,bound))
        U <= truerate && append!(points, t)
    end
    return points
end

thinning (generic function with 1 method)

Question 1¶

Simulation of a renewal process. Its intensity depends on the time elapsed since the last point, i.e. $\lambda_{t} = q(t-T_{N_{t-}})$, with $q(x) = K/(a+x)$. $K$ and $a$ are positive real parameters.

Here is the implementation of the corresponding functions intensity and intensitybound.

In [ ]:
function intensityQ1(t, history; K=1, a=1)
    isempty(history) ? lastpoint = 0 : lastpoint = last(history)
    s = t-lastpoint
    return K/(a+s)
end
function intensityboundQ1(t, history; K=1, a=1)
    isempty(history) ? lastpoint = 0 : lastpoint = last(history)
    s = t-lastpoint
    return (K/(a+s), Inf)
end

intensityboundQ1 (generic function with 1 method)

Here is how to call the thinning simulation of the model.

In [ ]:
K=2; a=2; tmax = 20;
points = thinning(tmax, (t,history)->intensityQ1(t, history; K=K, a=a), (t,history)->intensityboundQ1(t, history; K=K, a=a))

5-element Vector{Real}:
  8.515874360074603
  8.890720529968261
 10.143785017869178
 13.995783357566346
 14.707134526905278

Here are two functions to simplify notation.

In [ ]:
thinningQ1(tmax; K=1, a=1) = thinning(tmax, (t,history)->intensityQ1(t, history; K=K, a=a), (t,history)->intensityboundQ1(t, history; K=K, a=a))
graphQ1(tmax, points; K=1, a=1) = pointprocess_withintensity(0, tmax, points, (t,history)->intensityQ1(t, history; K=K, a=a))
In [ ]:
tmax=50.

K=1; a=1;
points = thinningQ1(tmax; K=K, a=a)
graphQ1(tmax, points; K=K, a=a)
display(title!("K="*string(K)*", a="*string(a)))
K=3; a=1;
points = thinningQ1(tmax; K=K, a=a)
graphQ1(tmax, points; K=K, a=a)
display(title!("K="*string(K)*", a="*string(a)))
K=1; a=.5;
points = thinningQ1(tmax; K=K, a=a)
graphQ1(tmax, points; K=K, a=a)
display(title!("K="*string(K)*", a="*string(a)))
K=1; a=2;
points = thinningQ1(tmax; K=K, a=a)
graphQ1(tmax, points; K=K, a=a)
display(title!("K="*string(K)*", a="*string(a)))

The parameter $K$ mainly controls the mean number of points whereas $a$ controls the clustering: $a$ small implies high clustering.

Question 2¶

Simulation of a linear Wold process. Its intensity depends on the time elapsed since the last point and the interval between the last two points, i.e. $\lambda_{t} = \mu+\alpha_{1}(t-T_{N_{t-}}) + \alpha_{2}(T_{N_{t-}} - T_{N_{t-}-1})$, with the convention that $T_{-1}=0$. $\mu$, $\alpha_{1}$ and $\alpha_{2}$ are positive real parameters.

Here is the implementation of the corresponding functions intensity and intensitybound.

In [ ]:
function intensityQ2(t, history; mu=1, alpha1=1, alpha2=1)
    if length(history)==0
        twolastpoints = [0,0]
    elseif length(history)==1
        twolastpoints = [0,last(history)] 
    else
        twolastpoints = last(history,2)
    end
    s1 = t-last(twolastpoints)
    s2 = diff(twolastpoints)[1]
    return mu + alpha1*s1 + alpha2*s2
end
function intensityboundQ2(t, history; mu=1, alpha1=1, alpha2=1)
    if length(history)==0
        twolastpoints = [0,0]
    elseif length(history)==1
        twolastpoints = [0,last(history)] 
    else
        twolastpoints = last(history,2)
    end
    tbound = t+1
    s1bound = tbound-last(twolastpoints)
    s2 = diff(twolastpoints)[1]
    return (mu + alpha1*s1bound + alpha2*s2, tbound)
end

intensityboundQ2 (generic function with 1 method)

Here is how to call the thinning simulation of the model.

In [ ]:
mu=1; alpha1=1; alpha2=1; tmax = 20;
points = thinning(tmax, (t,history)->intensityQ2(t, history; mu=mu, alpha1=alpha1, alpha2=alpha2), (t,history)->intensityboundQ2(t, history; mu=mu, alpha1=alpha1, alpha2=alpha2))

37-element Vector{Real}:
  0.553088990277999
  0.824341109663302
  1.9188777037418356
  2.45955736995026
  3.1933187933546083
  3.3681597398293897
  3.4713465701326145
  4.002077522103385
  4.936640016423601
  5.47147419008879
  ⋮
 17.065232383132095
 17.12923985667036
 17.50784394718668
 18.10855055028872
 18.549716856353715
 18.726905772781564
 18.877959688498436
 19.21935635688563
 19.39282671970031

Here are two functions to simplify notation.

In [ ]:
thinningQ2(tmax; mu=1, alpha1=1, alpha2=1) = thinning(tmax, (t,history)->intensityQ2(t, history; mu=mu, alpha1=alpha1, alpha2=alpha2), (t,history)->intensityboundQ2(t, history; mu=mu, alpha1=alpha1, alpha2=alpha2))
graphQ2(tmax, points; mu=1, alpha1=1, alpha2=1) = pointprocess_withintensity(0, tmax, points, (t,history)->intensityQ2(t, history; mu=mu, alpha1=alpha1, alpha2=alpha2))

graphQ2 (generic function with 1 method)
In [ ]:
tmax=20.

mu=1; alpha1=1; alpha2=1;
points = thinningQ2(tmax; mu=mu, alpha1=alpha1, alpha2=alpha2)
graphQ2(tmax, points; mu=mu, alpha1=alpha1, alpha2=alpha2)
display(title!("mu="*string(mu)*", α₁="*string(alpha1)*", α₂="*string(alpha2)))
mu=3; alpha1=1; alpha2=1;
points = thinningQ2(tmax; mu=mu, alpha1=alpha1, alpha2=alpha2)
graphQ2(tmax, points; mu=mu, alpha1=alpha1, alpha2=alpha2)
display(title!("mu="*string(mu)*", α₁="*string(alpha1)*", α₂="*string(alpha2)))
mu=1; alpha1=.5; alpha2=2;
points = thinningQ2(tmax; mu=mu, alpha1=alpha1, alpha2=alpha2)
graphQ2(tmax, points; mu=mu, alpha1=alpha1, alpha2=alpha2)
display(title!("mu="*string(mu)*", α₁="*string(alpha1)*", α₂="*string(alpha2)))
mu=1; alpha1=2; alpha2=.5;
points = thinningQ2(tmax; mu=mu, alpha1=alpha1, alpha2=alpha2)
graphQ2(tmax, points; mu=mu, alpha1=alpha1, alpha2=alpha2)
display(title!("mu="*string(mu)*", α₁="*string(alpha1)*", α₂="*string(alpha2)))
In [ ]: