diff --git a/docs/Project.toml b/docs/Project.toml
index 9ca9240..f8799cd 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -8,6 +8,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PointProcesses = "af0b7596-9bb0-472a-a012-63904f2b4c55"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/docs/examples/Basics.jl b/docs/examples/Basics.jl
index 11c7bbf..f954290 100644
--- a/docs/examples/Basics.jl
+++ b/docs/examples/Basics.jl
@@ -157,7 +157,7 @@ p_waiting
 # Let's fit the actual model now using PointProcesses.jl and run some formal statistical tests.
 
 # ### Fitting a Homogeneous Poisson Process
-pp_model = fit(PoissonProcess{Float64,Dirac{Nothing}}, h)
+pp_model = fit(PoissonProcess{Float64,NoMarks}, h)
 
 println("Estimated rate λ̂ = ", pp_model.λ) # hide
 
diff --git a/docs/examples/CustomProcesses.jl b/docs/examples/CustomProcesses.jl
new file mode 100644
index 0000000..c8119e5
--- /dev/null
+++ b/docs/examples/CustomProcesses.jl
@@ -0,0 +1,291 @@
+# # Defining Custom Processes
+
+# Although several processes are already implemented and more are to come, we might need
+# to define our own custom processes.  
+# Here we will explain the interface for defining new processes, starting with a simple case
+# and building up to more complex cases (which should not be too complicated).
+
+# ## Basic interface
+
+# For this example we will implement the two-state model from [Selva2022](@cite).
+# This model takes three parameters, λₗ, λₕ and τ. The process starts
+# as a Poisson model with intensity λ = λₗ. Whenever an event occurs, the process
+# immediately switches to a Poisson process with intensity λ = λₕ for τ
+# units of time. If an event occurs during this high intensity phase, the "counter" τ simply
+# resets. You can visualize it below.
+
+# ![Adapted from Selva (2022)](../assets/selva2022.png)
+
+# If the history ℋₜ up to time t contains the events t₁ < t₂ < ... < tₙ, then the equation
+# for the conditional intensity is  
+
+#     λ(t | ℋₜ) = λₗ + 𝟙(tₙ > t - τ) (λₕ - λₗ),
+
+# where 𝟙 is the indicator function. Notice that tₙ < t by definition of ℋₜ.
+
+# We will use a different parametrization, namely λ = λₗ and Δ = λₕ - λₗ, so the constraint
+# λₕ > λₗ becomes Δ > 0. 
+
+# The first step in our implementation is to define this new struct by subtyping
+# the `AbstractUnivariateProcess` type (we will get to multivariate processes).
+
+using PointProcesses
+using Distributions
+using Plots
+using StatsAPI
+using Optim
+
+import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess
+
+struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess
+    λ::R
+    Δ::R
+    τ::R
+    mark_dist::D
+end
+
+# > Note: Subtype `AbstractUnivariateProcess`, not just `AbstractPointProcess`, this way we can use `TwoStateModel` to build multivariate processes with `IndependentMultivariateProcess`.
+
+# `PointProcessMarkDistribution` encompasses `NoMarks` and any distribution in package
+# `Distributions.jl`. In the next chapter we will discuss custom distributions as well.
+
+# > Even if you plan on only using a non-marked process, you should keep the `mark_dist` field. If you really want you can define an inner constructor like `TwoStateModel(λ, Δ, τ) → new(λ, Δ, τ, Nomarks())`.
+
+# Right now you cannot do much more than accessing the fields. The only methods you have
+# defined are `ndims`, `DensityKind` and `mark_distribution`. Not too interesting.  
+
+# We will begin by implementing the methods `ground_intensity`, `integrated_ground_intensity` and
+# `ground_intensity_bound`. They are sufficient for almost all methods to function properly.
+
+function PointProcesses.ground_intensity(pp::TwoStateModel, t, h::History)
+    if isempty(h) || t <= h.times[1]
+        return pp.λ
+    else
+        t_n = h.times[searchsortedfirst(h.times, t) - 1] # Last event time before `t`
+        in_high_phase = t < t_n + pp.τ
+    end
+    return pp.λ + (in_high_phase * pp.Δ)
+end
+
+function PointProcesses.ground_intensity_bound(pp::TwoStateModel, t, h::History)
+    return (pp.λ + pp.Δ, Inf) # Also not the most efficient implementation, but works
+end
+
+function PointProcesses.integrated_ground_intensity(pp::TwoStateModel, h, a, b)
+    integral = pp.λ * (b - a)
+    interval_events = event_times(h, a, b)
+    if length(interval_events) > 0
+        for i in 1:(length(interval_events) - 1)
+            integral += pp.Δ * min(interval_events[i + 1] - interval_events[i], pp.τ)
+        end
+        integral += pp.Δ * min(b - interval_events[end], pp.τ)
+    end
+    return integral
+end
+
+# > Notice that in `ground_intensity` we searched for the last event that occurred before time `t`. We need to assume that `h` might be the full history, and not only the history up to time `t`.
+
+# With this we can now call many other methods: `intensity`, `log_intensity`, `logdensityof`, `time_change`
+# and `simulate`. Here are the dependencies:
+# - `intensity` → `ground_intensity` + `densityof` from `mark_dist`
+# - `log_intensity` → `intensity`
+# - `logdensityof` → `integrated_ground_intensity` + `logdensityof`
+# - `time_change` → `integrated_ground_intensity`
+# - `simulate` → `ground_intensity` + `ground_intensity_bound`
+
+pp = TwoStateModel(1.0, 2.0, 0.5, Normal())
+h = simulate(pp, 0.0, 20.0)
+
+println("intensity(pp, m, t, h) = $(intensity(pp, 0.0, 0.0, h))")
+println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))")
+println("logdensityof(pp, h) = $(logdensityof(pp, h))")
+println("time_change(pp, h) = $(time_change(h, pp))")
+println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))")
+
+# Let's plot our simulated process `h` along with its ground intensity.
+
+sim_plot = plot(
+    event_times(h),
+    fill(1.0, nb_events(h));
+    line=:stem,
+    grid=false,
+    label=false,
+    xlabel="Time",
+    ylabel="Events",
+    yaxis=false,
+    xlim=(min_time(h), max_time(h)),
+)
+
+xs = LinRange(min_time(h), max_time(h), 1000)
+intensity_plot = plot(
+    xs,
+    [ground_intensity(pp, x, h) for x in xs],
+    ;
+    xaxis=false,
+    xgrid=false,
+    label=false,
+    ylabel="Ground Intensity",
+    xlim=(min_time(h), max_time(h)),
+)
+
+plot(
+    intensity_plot,
+    sim_plot;
+    title=["Two State Model Simulation" ""],
+    layout=grid(2, 1; heights=[0.7, 0.3]),
+)
+
+# Lastly, we can define the `fit` method for estimating the parameters from an observed
+# history. `logdensityof` is simply the log-likelihood of our process, so for a quick
+# and simple example we keep the mark distribution fixed and only optimize the core
+# process parameters.
+
+function StatsAPI.fit(
+    ::Type{TwoStateModel{R1,NoMarks}}, h::History, init_params::Vector{R2}
+) where {R1<:Real,R2<:Real}
+    objective(params) = -logdensityof(TwoStateModel(params..., NoMarks()), h)
+
+    lower_bound = [0.0, 0.0, 0.0]
+    upper_bound = [Inf, Inf, Inf]
+    result = optimize(objective, lower_bound, upper_bound, init_params)
+    optimal = Optim.minimizer(result)
+
+    return TwoStateModel(optimal..., NoMarks())
+end
+
+# > It is tricky to properly estimate the parameters from this process. This implementation is conceptually correct, but it should be improved in real applications.
+
+# With the `fit`, `simulate` and `time_change` methods, we can perform
+# goodness-of-fit tests for our model.
+
+function random_params()
+    λ = 10.0 + rand() * 10
+    Δ = λ * rand()
+    τ = inv(λ + Δ) * (0.5 + rand())
+    return [λ, Δ, τ]
+end
+
+true_params = random_params()
+h_unknown = simulate(TwoStateModel(true_params..., NoMarks()), 0.0, 100.0)
+
+init_params = random_params()
+pp_estimated = fit(TwoStateModel{Float64,NoMarks}, h_unknown, init_params)
+estimated_params = [pp_estimated.λ, pp_estimated.Δ, pp_estimated.τ] # hide
+
+println("True and estimated parameters:") # hide
+println("λ: $(true_params[1]) - $(estimated_params[1])") #hide
+println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide
+println("τ : $(true_params[3]) - $(estimated_params[3])") # hide
+
+pp_to_test = TwoStateModel(20.0, 30.0, 1.0 / 30.0, NoMarks())
+# n_sims=100 for the sake of this example. In real applications, this should be higher
+test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown; n_sims=100)
+
+println("p-value for hypothesis test: $(pvalue(test_result))") # hide
+
+# > The statistics for goodness of fit tests currently implemented are suited only for testing the distribution of event times. They will run for marked processes, but the test results are only related to the event times.
+
+# ## Custom Mark Distributions
+
+# Mark distributions can also be customized. Analogously to processes, we need
+# to subtype the `AbstractMarkDistribution` type. Lets make a mark distribution
+# that changes with times.
+
+struct LinearTimeNormal{R<:Real} <: AbstractMarkDistribution
+    L::R
+    H::R
+end
+
+# There are five methods that can be used with an `AbstractMarkDistribution`: `mark_distribution`,
+# `sample_mark`, `DensityInterface.densityof`, `Base.eltype` and `StatsAPI.fit`. The number of
+# methods we need to implement will depend on what `mark_distribution` returns.
+
+# `mark_distribution` is used to return the distribution of the marks at a specific time `t` after
+# some history `h`. Say `mark_distribution(md::AbstractMarkDistribution, t, h::History)` is of
+# type `T`.
+# - `rand(::AbstractRNG, ::T)` implemented ⇒ `sample_mark(md, t, h)` works
+# - `eltype(::T)` implemented ⇒ `eltype(md)` works
+# - `densityof(::T, m)` implemented ⇒ `densityof(md, t, h, m)` works
+# - `fit` must always be implemented
+
+function PointProcesses.mark_distribution(md::LinearTimeNormal, t, h::History)
+    time_ratio = (t - min_time(h)) / duration(h)
+    return Normal(md.L + (md.H - md.L) * time_ratio, 1)
+end
+
+# > It is important to think of the case where the history `h` is empty. If your custom distribution depends on past events, there should be a specific branch covering the case `isempty(h)`. 
+
+# In our case `mark_distribution` returns a `Distribution`, so `sample_mark`, `eltype` and
+# `densityof` are already take care of. To have all the functionality from the standard
+# mark distributions, we would only need to implement `fit`.
+
+pp = TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0))
+
+println(
+    "TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) = $(TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)))",
+)
+println(
+    "PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)) = $(PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)))",
+)
+println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))")
+println("logdensityof(pp, h) = $(logdensityof(pp, h))")
+println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))")
+
+# ## Multivariate Processes
+
+# Lastly we can define a multivariate process composed of multiple independent
+# univariate processes. If a method is available for all individual processes,
+# it will be also available for the multivariate processes.  
+# All methods now return a vector, one for each dimension (each univariate process),
+# and we can add an integer argument `d` for accessing the value for a specific
+# dimension.  
+
+imp = IndependentMultivariateProcess([
+    TwoStateModel(1.2, 2.0, 0.3, Normal()),
+    TwoStateModel(0.8, 2.0, 0.4, NoMarks()),
+    PoissonProcess(2.0), # Any type of process can be added
+])
+
+sim = simulate(imp, 0.0, 10.0)
+
+println("ground_intensity(imp, 0.0, sim) → $(ground_intensity(imp, 0.0, sim))")
+println("ground_intensity(imp, 0.0, sim, 1) → $(ground_intensity(imp, 0.0, sim, 1))")
+println("intensity(imp, 0.0, 0.0, sim) → $(intensity(imp, 0.0, 0.0, sim))")
+println("intensity(imp, 0.0, nothing, sim, 3) → $(intensity(imp, nothing, 0.0, sim, 3))")
+
+# Lets end this tutorial with a nice plot of this process. 
+
+λ(t) = ground_intensity(imp, t, sim) # Returns a vector with the 3 individual intensities
+
+events_plot = scatter()
+for d in 1:ndims(sim)
+    scatter!(
+        events_plot,
+        event_times(sim, d),
+        fill(d, nb_events(sim, d));
+        label=nothing,
+        markersize=3,
+    )
+end
+plot!(
+    events_plot;
+    yaxis=false,
+    xlabel="Time",
+    xlim=(min_time(sim), max_time(sim)),
+    ylim=(0.0, 3.2),
+)
+
+xs = LinRange(min_time(sim), max_time(sim), 1000)
+intensities_plot = plot()
+for d in 1:ndims(sim)
+    plot!(intensities_plot, xs, getindex.(λ.(xs), d); label="Dimension $d")
+end
+plot!(intensities_plot, xs, sum.(λ.(xs)); label="Total ground intensity")
+plot!(intensities_plot; xaxis=false, xlim=(min_time(sim), max_time(sim)), legend=:topleft)
+
+plot(
+    intensities_plot,
+    events_plot;
+    title="Event times and ground intensities",
+    layout=grid(2, 1; heights=[0.7, 0.3]),
+)
diff --git a/docs/examples/Hawkes.jl b/docs/examples/Hawkes.jl
index 44d37b5..59a810f 100644
--- a/docs/examples/Hawkes.jl
+++ b/docs/examples/Hawkes.jl
@@ -106,9 +106,7 @@ tod = sort(tod_raw .+ (1:length(tod_raw)) .* 1e-6)
 h_day = History(tod, 0.0, 1440.0) # build a "history" on [0, 1440] minutes
 nbins = 96  # 96 bins = 15-minute bins
 pp_day = fit(
-    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
-    h_day,
-    nbins,
+    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h_day, nbins
 )
 
 λ_avg(u) = pp_day.intensity_function(u) / n_days
diff --git a/docs/examples/Inhomogeneous.jl b/docs/examples/Inhomogeneous.jl
index 7e0bafa..6df7238 100644
--- a/docs/examples/Inhomogeneous.jl
+++ b/docs/examples/Inhomogeneous.jl
@@ -70,7 +70,7 @@ scatter!(
 # The simplest approach is to bin the data and estimate a constant rate in each bin:
 
 pp_piecewise = fit(
-    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
+    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks},
     h,
     20,  # number of bins
 )
diff --git a/docs/src/api.md b/docs/src/api.md
index 27d1f26..ad2454c 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -38,6 +38,12 @@ time_change
 split_into_chunks
 ```
 
+
+### Interface
+
+```@docs
+```
+
 ## Point processes
 
 ```@docs
@@ -47,20 +53,24 @@ AbstractMultivariateProcess
 BoundedPointProcess
 IndependentMultivariateProcess
 Base.ndims(::AbstractPointProcess)
+eltype
 ```
 
-### Intensity
+### Mark Distributions
 
 ```@docs
-intensity
-ground_intensity
-log_intensity
+AbstractMarkDistribution
+NoMarks
 ```
 
-### Marks
+### Intensity
 
 ```@docs
+intensity
+ground_intensity
+log_intensity
 mark_distribution
+densityof
 ```
 
 ### Simulation
@@ -68,6 +78,7 @@ mark_distribution
 ```@docs
 simulate_ogata
 simulate
+sample_mark
 ```
 
 ### Inference
diff --git a/docs/src/assets/selva2022.png b/docs/src/assets/selva2022.png
new file mode 100644
index 0000000..7d6391c
Binary files /dev/null and b/docs/src/assets/selva2022.png differ
diff --git a/docs/src/examples/Basics.md b/docs/src/examples/Basics.md
index 13a5eb8..c565a5a 100644
--- a/docs/src/examples/Basics.md
+++ b/docs/src/examples/Basics.md
@@ -98,7 +98,7 @@ counts to be roughly Poisson distributed with rate λ * Δt in each bin of width
 
 ````@example Basics
 bin_width = 1.0
-bins = collect(h.tmin:bin_width:h.tmax)           # bin edges
+bins = collect((h.tmin):bin_width:(h.tmax))           # bin edges
 counts = fit(Histogram, h.times, bins).weights;    # counts per bin
 nothing #hide
 ````
@@ -185,7 +185,7 @@ Let's fit the actual model now using PointProcesses.jl and run some formal stati
 ### Fitting a Homogeneous Poisson Process
 
 ````@example Basics
-pp_model = fit(PoissonProcess{Float64,Dirac{Nothing}}, h)
+pp_model = fit(PoissonProcess{Float64,NoMarks}, h)
 
 println("Estimated rate λ̂ = ", pp_model.λ) # hide
 ````
diff --git a/docs/src/examples/CustomProcesses.md b/docs/src/examples/CustomProcesses.md
new file mode 100644
index 0000000..8104634
--- /dev/null
+++ b/docs/src/examples/CustomProcesses.md
@@ -0,0 +1,322 @@
+```@meta
+EditURL = "../../examples/CustomProcesses.jl"
+```
+
+# Defining Custom Processes
+
+Although several processes are already implemented and more are to come, we might need
+to define our own custom processes.
+Here we will explain the interface for defining new processes, starting with a simple case
+and building up to more complex cases (which should not be too complicated).
+
+## Basic interface
+
+For this example we will implement the two-state model from [Selva2022](@cite).
+This model takes three parameters, λₗ, λₕ and τ. The process starts
+as a Poisson model with intensity λ = λₗ. Whenever an event occurs, the process
+immediately switches to a Poisson process with intensity λ = λₕ for τ
+units of time. If an event occurs during this high intensity phase, the "counter" τ simply
+resets. You can visualize it below.
+
+![Adapted from Selva (2022)](../assets/selva2022.png)
+
+If the history ℋₜ up to time t contains the events t₁ < t₂ < ... < tₙ, then the equation
+for the conditional intensity is
+
+    λ(t | ℋₜ) = λₗ + 𝟙(tₙ > t - τ) (λₕ - λₗ),
+
+where 𝟙 is the indicator function. Notice that tₙ < t by definition of ℋₜ.
+
+We will use a different parametrization, namely λ = λₗ and Δ = λₕ - λₗ, so the constraint
+λₕ > λₗ becomes Δ > 0.
+
+The first step in our implementation is to define this new struct by subtyping
+the `AbstractUnivariateProcess` type (we will get to multivariate processes).
+
+````@example CustomProcesses
+using PointProcesses
+using Distributions
+using Plots
+using StatsAPI
+using Optim
+
+import PointProcesses: NoMarks, AbstractMarkDistribution, AbstractUnivariateProcess
+
+struct TwoStateModel{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess
+    λ::R
+    Δ::R
+    τ::R
+    mark_dist::D
+end
+````
+
+> Note: Subtype `AbstractUnivariateProcess`, not just `AbstractPointProcess`, this way we can use `TwoStateModel` to build multivariate processes with `IndependentMultivariateProcess`.
+
+`PointProcessMarkDistribution` encompasses `NoMarks` and any distribution in package
+`Distributions.jl`. In the next chapter we will discuss custom distributions as well.
+
+> Even if you plan on only using a non-marked process, you should keep the `mark_dist` field. If you really want you can define an inner constructor like `TwoStateModel(λ, Δ, τ) → new(λ, Δ, τ, Nomarks())`.
+
+Right now you cannot do much more than accessing the fields. The only methods you have
+defined are `ndims`, `DensityKind` and `mark_distribution`. Not too interesting.
+
+We will begin by implementing the methods `ground_intensity`, `integrated_ground_intensity` and
+`ground_intensity_bound`. They are sufficient for almost all methods to function properly.
+
+````@example CustomProcesses
+function PointProcesses.ground_intensity(pp::TwoStateModel, t, h::History)
+    if isempty(h) || t <= h.times[1]
+        return pp.λ
+    else
+        t_n = h.times[searchsortedfirst(h.times, t) - 1] # Last event time before `t`
+        in_high_phase = t < t_n + pp.τ
+    end
+    return pp.λ + (in_high_phase * pp.Δ)
+end
+
+function PointProcesses.ground_intensity_bound(pp::TwoStateModel, t, h::History)
+    return (pp.λ + pp.Δ, Inf) # Also not the most efficient implementation, but works
+end
+
+function PointProcesses.integrated_ground_intensity(pp::TwoStateModel, h, a, b)
+    integral = pp.λ * (b - a)
+    interval_events = event_times(h, a, b)
+    if length(interval_events) > 0
+        for i in 1:(length(interval_events) - 1)
+            integral += pp.Δ * min(interval_events[i + 1] - interval_events[i], pp.τ)
+        end
+        integral += pp.Δ * min(b - interval_events[end], pp.τ)
+    end
+    return integral
+end
+````
+
+> Notice that in `ground_intensity` we searched for the last event that occurred before time `t`. We need to assume that `h` might be the full history, and not only the history up to time `t`.
+
+With this we can now call many other methods: `intensity`, `log_intensity`, `logdensityof`, `time_change`
+and `simulate`. Here are the dependencies:
+- `intensity` → `ground_intensity` + `densityof` from `mark_dist`
+- `log_intensity` → `intensity`
+- `logdensityof` → `integrated_ground_intensity` + `logdensityof`
+- `time_change` → `integrated_ground_intensity`
+- `simulate` → `ground_intensity` + `ground_intensity_bound`
+
+````@example CustomProcesses
+pp = TwoStateModel(1.0, 2.0, 0.5, Normal())
+h = simulate(pp, 0.0, 20.0)
+
+println("intensity(pp, m, t, h) = $(intensity(pp, 0.0, 0.0, h))")
+println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))")
+println("logdensityof(pp, h) = $(logdensityof(pp, h))")
+println("time_change(pp, h) = $(time_change(h, pp))")
+println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))")
+````
+
+Let's plot our simulated process `h` along with its ground intensity.
+
+````@example CustomProcesses
+sim_plot = plot(
+    event_times(h),
+    fill(1.0, nb_events(h));
+    line=:stem,
+    grid=false,
+    label=false,
+    xlabel="Time",
+    ylabel="Events",
+    yaxis=false,
+    xlim=(min_time(h), max_time(h)),
+)
+
+xs = LinRange(min_time(h), max_time(h), 1000)
+intensity_plot = plot(
+    xs,
+    [ground_intensity(pp, x, h) for x in xs],
+    ;
+    xaxis=false,
+    xgrid=false,
+    label=false,
+    ylabel="Ground Intensity",
+    xlim=(min_time(h), max_time(h)),
+)
+
+plot(
+    intensity_plot,
+    sim_plot;
+    title=["Two State Model Simulation" ""],
+    layout=grid(2, 1; heights=[0.7, 0.3]),
+)
+````
+
+Lastly, we can define the `fit` method for estimating the parameters from an observed
+history. `logdensityof` is simply the log-likelihood of our process, so for a quick
+and simple example we keep the mark distribution fixed and only optimize the core
+process parameters.
+
+````@example CustomProcesses
+function StatsAPI.fit(
+    ::Type{TwoStateModel{R1,NoMarks}}, h::History, init_params::Vector{R2}
+) where {R1<:Real,R2<:Real}
+    objective(params) = -logdensityof(TwoStateModel(params..., NoMarks()), h)
+
+    lower_bound = [0.0, 0.0, 0.0]
+    upper_bound = [Inf, Inf, Inf]
+    result = optimize(objective, lower_bound, upper_bound, init_params)
+    optimal = Optim.minimizer(result)
+
+    return TwoStateModel(optimal..., NoMarks())
+end
+````
+
+> It is tricky to properly estimate the parameters from this process. This implementation is conceptually correct, but it should be improved in real applications.
+
+With the `fit`, `simulate` and `time_change` methods, we can perform
+goodness-of-fit tests for our model.
+
+````@example CustomProcesses
+function random_params()
+    λ = 10.0 + rand() * 10
+    Δ = λ * rand()
+    τ = inv(λ + Δ) * (0.5 + rand())
+    return [λ, Δ, τ]
+end
+
+true_params = random_params()
+h_unknown = simulate(TwoStateModel(true_params..., NoMarks()), 0.0, 100.0)
+
+init_params = random_params()
+pp_estimated = fit(TwoStateModel{Float64,NoMarks}, h_unknown, init_params)
+estimated_params = [pp_estimated.λ, pp_estimated.Δ, pp_estimated.τ] # hide
+
+println("True and estimated parameters:") # hide
+println("λ: $(true_params[1]) - $(estimated_params[1])") #hide
+println("Δ: $(true_params[2]) - $(estimated_params[2])") #hide
+println("τ : $(true_params[3]) - $(estimated_params[3])") # hide
+
+pp_to_test = TwoStateModel(20.0, 30.0, 1.0 / 30.0, NoMarks())
+````
+
+n_sims=100 for the sake of this example. In real applications, this should be higher
+
+````@example CustomProcesses
+test_result = MonteCarloTest(KSDistance{Exponential}, pp_to_test, h_unknown; n_sims=100)
+
+println("p-value for hypothesis test: $(pvalue(test_result))") # hide
+````
+
+> The statistics for goodness of fit tests currently implemented are suited only for testing the distribution of event times. They will run for marked processes, but the test results are only related to the event times.
+
+## Custom Mark Distributions
+
+Mark distributions can also be customized. Analogously to processes, we need
+to subtype the `AbstractMarkDistribution` type. Lets make a mark distribution
+that changes with times.
+
+````@example CustomProcesses
+struct LinearTimeNormal{R<:Real} <: AbstractMarkDistribution
+    L::R
+    H::R
+end
+````
+
+There are five methods that can be used with an `AbstractMarkDistribution`: `mark_distribution`,
+`sample_mark`, `DensityInterface.densityof`, `Base.eltype` and `StatsAPI.fit`. The number of
+methods we need to implement will depend on what `mark_distribution` returns.
+
+`mark_distribution` is used to return the distribution of the marks at a specific time `t` after
+some history `h`. Say `mark_distribution(md::AbstractMarkDistribution, t, h::History)` is of
+type `T`.
+- `rand(::AbstractRNG, ::T)` implemented ⇒ `sample_mark(md, t, h)` works
+- `eltype(::T)` implemented ⇒ `eltype(md)` works
+- `densityof(::T, m)` implemented ⇒ `densityof(md, t, h, m)` works
+- `fit` must always be implemented
+
+````@example CustomProcesses
+function PointProcesses.mark_distribution(md::LinearTimeNormal, t, h::History)
+    time_ratio = (t - min_time(h)) / duration(h)
+    return Normal(md.L + (md.H - md.L) * time_ratio, 1)
+end
+````
+
+> It is important to think of the case where the history `h` is empty. If your custom distribution depends on past events, there should be a specific branch covering the case `isempty(h)`.
+
+In our case `mark_distribution` returns a `Distribution`, so `sample_mark`, `eltype` and
+`densityof` are already take care of. To have all the functionality from the standard
+mark distributions, we would only need to implement `fit`.
+
+````@example CustomProcesses
+pp = TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0))
+
+println(
+    "TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)) = $(TwoStateModel(1.0, 2.0, 3.0, LinearTimeNormal(2.0, 4.0)))",
+)
+println(
+    "PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)) = $(PoissonProcess(2.0, LinearTimeNormal(2.0, 4.0)))",
+)
+println("log_intensity(pp, m, t, h) = $(log_intensity(pp, 0.0, 0.0, h))")
+println("logdensityof(pp, h) = $(logdensityof(pp, h))")
+println("simulate(pp, 0.0, 10.0) = $(simulate(pp, 0.0, 10.0))")
+````
+
+## Multivariate Processes
+
+Lastly we can define a multivariate process composed of multiple independent
+univariate processes. If a method is available for all individual processes,
+it will be also available for the multivariate processes.
+All methods now return a vector, one for each dimension (each univariate process),
+and we can add an integer argument `d` for accessing the value for a specific
+dimension.
+
+````@example CustomProcesses
+imp = IndependentMultivariateProcess([
+    TwoStateModel(1.2, 2.0, 0.3, Normal()),
+    TwoStateModel(0.8, 2.0, 0.4, NoMarks()),
+    PoissonProcess(2.0), # Any type of process can be added
+])
+
+sim = simulate(imp, 0.0, 10.0)
+
+println("ground_intensity(imp, 0.0, sim) → $(ground_intensity(imp, 0.0, sim))")
+println("ground_intensity(imp, 0.0, sim, 1) → $(ground_intensity(imp, 0.0, sim, 1))")
+println("intensity(imp, 0.0, 0.0, sim) → $(intensity(imp, 0.0, 0.0, sim))")
+println("intensity(imp, 0.0, nothing, sim, 3) → $(intensity(imp, nothing, 0.0, sim, 3))")
+````
+
+Lets end this tutorial with a nice plot of this process.
+
+````@example CustomProcesses
+λ(t) = ground_intensity(imp, t, sim) # Returns a vector with the 3 individual intensities
+
+events_plot = scatter()
+for d in 1:ndims(sim)
+    scatter!(
+        events_plot,
+        event_times(sim, d),
+        fill(d, nb_events(sim, d));
+        label=nothing,
+        markersize=3,
+    )
+end
+plot!(
+    events_plot;
+    yaxis=false,
+    xlabel="Time",
+    xlim=(min_time(sim), max_time(sim)),
+    ylim=(0.0, 3.2),
+)
+
+xs = LinRange(min_time(sim), max_time(sim), 1000)
+intensities_plot = plot()
+for d in 1:ndims(sim)
+    plot!(intensities_plot, xs, getindex.(λ.(xs), d); label="Dimension $d")
+end
+plot!(intensities_plot, xs, sum.(λ.(xs)); label="Total ground intensity")
+plot!(intensities_plot; xaxis=false, xlim=(min_time(sim), max_time(sim)), legend=:topleft)
+
+plot(
+    intensities_plot,
+    events_plot;
+    title="Event times and ground intensities",
+    layout=grid(2, 1; heights=[0.7, 0.3]),
+)
+````
+
diff --git a/docs/src/examples/Hawkes.md b/docs/src/examples/Hawkes.md
index 6af2ed0..05c7237 100644
--- a/docs/src/examples/Hawkes.md
+++ b/docs/src/examples/Hawkes.md
@@ -59,7 +59,7 @@ dt = DateTime("12/22 3:09 pm", fmt)
 
 function parse_timestamp_min(s::AbstractString; year=2024)
     dt = DateTime(s, dateformat"m/d I:M p")
-    DateTime(year, month(dt), day(dt), hour(dt), minute(dt))
+    return DateTime(year, month(dt), day(dt), hour(dt), minute(dt))
 end
 
 data.TimestampDT = parse_timestamp_min.(String.(data.Timestamp))
@@ -89,7 +89,7 @@ function eventplot(
     xlabel="Time (minutes)",
     ylabel="Events",
 )
-    scatter(
+    return scatter(
         event_times,
         ones(length(event_times));
         markershape=:vline,
@@ -135,9 +135,7 @@ tod = sort(tod_raw .+ (1:length(tod_raw)) .* 1e-6)
 h_day = History(tod, 0.0, 1440.0) # build a "history" on [0, 1440] minutes
 nbins = 96  # 96 bins = 15-minute bins
 pp_day = fit(
-    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
-    h_day,
-    nbins,
+    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks}, h_day, nbins
 )
 
 λ_avg(u) = pp_day.intensity_function(u) / n_days
@@ -180,8 +178,9 @@ Finally, we can visualize the fitted intensity function over time.
 ts = sort(data.t)
 
 function λ_hawkes(t::Real)
-    hawkes_model.μ +
-    sum((hawkes_model.α * exp(-hawkes_model.ω * (t - ti)) for ti in ts if ti < t); init=0.0)
+    return hawkes_model.μ + sum(
+        (hawkes_model.α * exp(-hawkes_model.ω * (t - ti)) for ti in ts if ti < t); init=0.0
+    )
 end
 
 u = range(0.0, maximum(ts) + 1.0; length=2000)
diff --git a/docs/src/examples/Inhomogeneous.md b/docs/src/examples/Inhomogeneous.md
index 901d21b..c625c37 100644
--- a/docs/src/examples/Inhomogeneous.md
+++ b/docs/src/examples/Inhomogeneous.md
@@ -92,7 +92,7 @@ The simplest approach is to bin the data and estimate a constant rate in each bi
 
 ````@example Inhomogeneous
 pp_piecewise = fit(
-    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
+    InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks},
     h,
     20,  # number of bins
 )
@@ -307,7 +307,7 @@ models = [
 ]
 
 println("\nModel Comparison (Negative Log-Likelihood):") # hide
-println("-" ^ 50) # hide
+println("-"^50) # hide
 for (name, model) in models # hide
     nll = compute_nll(model, h) # hide
     println("  $name: ", round(nll; digits=2)) # hide
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
index 8000823..9e5b630 100644
--- a/docs/src/refs.bib
+++ b/docs/src/refs.bib
@@ -110,4 +110,16 @@ @Book{Cox2018
   isbn      = {9780203743034},
   month     = dec,
   doi       = {10.1201/9780203743034},
+}
+
+@Article{Selva2022,
+  author    = {Jacopo Selva and Laura Sandri and Matteo Taroni and Roberto Sulpizio and Pablo Tierz and Antonio Costa},
+  journal   = {Science Advances},
+  title     = {A simple two-state model interprets temporal modulations in eruptive activity and enhances multivolcano hazard quantification},
+  year      = {2022},
+  month     = {nov},
+  number    = {44},
+  volume    = {8},
+  doi       = {10.1126/sciadv.abq4415},
+  publisher = {American Association for the Advancement of Science ({AAAS})},
 }
\ No newline at end of file
diff --git a/src/PointProcesses.jl b/src/PointProcesses.jl
index e37f5da..9315f5b 100644
--- a/src/PointProcesses.jl
+++ b/src/PointProcesses.jl
@@ -8,7 +8,8 @@ module PointProcesses
 # Imports
 
 using DensityInterface: DensityInterface, HasDensity, densityof, logdensityof
-using Distributions: Distributions, UnivariateDistribution, MultivariateDistribution
+using Distributions:
+    Distributions, Distribution, UnivariateDistribution, MultivariateDistribution
 using Distributions: Categorical, Exponential, Poisson, Uniform, Dirac, Gamma
 using Distributions: fit, suffstats, probs
 using Integrals: Integrals, IntegralProblem, solve, QuadGKJL
@@ -38,6 +39,10 @@ export min_time, max_time, min_mark, max_mark
 export nb_events, has_events, duration, ndims
 export time_change, split_into_chunks
 
+## Mark Distributions
+export AbstractMarkDistribution, PointProcessMarkDistribution, NoMarks
+export sample_mark
+
 ## Point processes
 
 export AbstractPointProcess, AbstractUnivariateProcess, AbstractMultivariateProcess
@@ -79,6 +84,7 @@ export PointProcessTest, BootstrapTest, MonteCarloTest
 
 ## General
 include("history.jl")
+include("mark_distributions.jl")
 include("abstract_point_process.jl")
 include("simulation.jl")
 include("bounded_point_process.jl")
diff --git a/src/abstract_point_process.jl b/src/abstract_point_process.jl
index 787b7b5..72810f2 100644
--- a/src/abstract_point_process.jl
+++ b/src/abstract_point_process.jl
@@ -55,7 +55,7 @@ For multivariate processes, it returns a vector of mark distributions for each d
 
 mark_distribution(pp, h, t, d) computes the mark distribution for a multivariate process `pp` at dimension `d`.
 """
-function mark_distribution end
+mark_distribution(pp::AbstractPointProcess, t, h) = mark_distribution(pp.mark_dist, t, h)
 
 """
     intensity(pp, m, t, h)
@@ -67,7 +67,9 @@ intensity(pp, h, t, d) computes the intensity for a multivariate process `pp` at
 
 The conditional intensity function `λ(t,m|h)` quantifies the instantaneous risk of an event with mark `m` occurring at time `t` after history `h`.
 """
-function intensity end
+function intensity(pp::AbstractPointProcess, m, t, h)
+    return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m)
+end
 
 """
     log_intensity(pp, m, t, h)
@@ -78,7 +80,7 @@ For multivariate processes, it returns a vector of log intensities for each dime
 log_intensity(pp, h, t, d) computes the log intensity for a multivariate process `pp` at dimension `d`.
 """
 function log_intensity(pp::AbstractPointProcess, m, t, h)
-    return log(ground_intensity(pp, t, h)) + logdensityof(mark_distribution(pp, t, h), m)
+    return log(intensity(pp, m, t, h))
 end
 
 ## Simulation
@@ -147,3 +149,20 @@ Fit a point process of type `PP` to one or several histories using maximum a pos
 Not implemented by default.
 """
 function fit_map end
+
+function time_change(h::History, pp::AbstractPointProcess)
+    Λ(t) = integrated_ground_intensity(pp, h, min_time(h), t)
+    return time_change(h, Λ)
+end
+
+"""
+    simulate([rng,] pp, tmin, tmax)
+
+Generate an event record from the point process `pp` on the interval [tmin, tmax).
+
+# Technical Remark
+Simulating events using types other than `Float64` may cause underflow issues, causing an error to be raised because of repeated event times. This is because the `rand` method from `Distributions.jl` always returns `Float64` values.
+"""
+function simulate(rng::AbstractRNG, pp::AbstractPointProcess, tmin, tmax)
+    return simulate_ogata(rng, pp, tmin, tmax)
+end
diff --git a/src/history.jl b/src/history.jl
index 68920f7..aa830bf 100644
--- a/src/history.jl
+++ b/src/history.jl
@@ -117,8 +117,12 @@ struct History{T<:Real,M,D}
 end
 
 function History(
-    times::Vector{Vector{R1}}, tmin::R2, tmax::R3, marks::Vector{Vector{M}}; check_args=true
-) where {R1<:Real,R2<:Real,R3<:Real,M}
+    times::AbstractVector{<:AbstractVector{R1}},
+    tmin::R2,
+    tmax::R3,
+    marks::AbstractVector{<:AbstractVector};
+    check_args=true,
+) where {R1<:Real,R2<:Real,R3<:Real}
     vec_times = vcat(times...)
     vec_marks = vcat(marks...)
     perm = sortperm(vec_times)
@@ -140,7 +144,7 @@ function History(
 end
 
 function History(
-    times::Vector{Vector{R1}}, tmin::R2, tmax::R3; check_args=true
+    times::AbstractVector{<:AbstractVector{R1}}, tmin::R2, tmax::R3; check_args=true
 ) where {R1<:Real,R2<:Real,R3<:Real}
     nots = [fill(nothing, length(times[i])) for i in 1:length(times)]
     return History(times, tmin, tmax, nots; check_args=check_args)
@@ -158,6 +162,17 @@ function History(; times, tmin, tmax, marks=nothing, dims=nothing, check_args=tr
     end
 end
 
+function History(tmin::R1, tmax::R2, N::Int=1) where {R1<:Real,R2<:Real}
+    R = promote_type(R1, R2)
+    return History(R[], tmin, tmax, [], [], N)
+end
+
+function History(h::History, d::Int)
+    times = event_times(h, d)
+    marks = event_marks(h, d)
+    return History(times, min_time(h), max_time(h), marks)
+end
+
 function Base.show(io::IO, h::History{T,M}) where {T,M}
     return print(
         io, "History{$T,$M} with $(nb_events(h)) events on interval [$(h.tmin), $(h.tmax))"
@@ -176,8 +191,8 @@ event_times(h::History) = h.times
 
 Return the sorted vector of event times for `h` in dimension `d`.
 """
-function event_times(h::History, d::Int)
-    return h.N == 1 && d == 1 ? h.times : (@view h.times[h.dims .== d])
+function event_times(h::History, d::Union{Int,Nothing})
+    return h.N == 1 ? h.times : (@view h.times[h.dims .== d])
 end
 
 """
@@ -215,7 +230,7 @@ event_marks(h::History) = h.marks
 
 Return the vector of event marks in dimension `d` of `h`, sorted according to their event times.
 """
-function event_marks(h::History, d::Int)
+function event_marks(h::History, d::Union{Int,Nothing})
     return h.N == 1 && d == 1 ? h.marks : (@view h.marks[h.dims .== d])
 end
 
@@ -354,17 +369,37 @@ duration(h::History) = max_time(h) - min_time(h)
 
 Add event `(t, m)` inside the interval `[h.tmin, h.tmax)` at the end of history `h`.
 
-With `check_args=true` (the default), the event must satisfy
-`h.tmin <= t < h.tmax`, occur at or after the last existing event time, and
-(if `h` is multivariate) lie in dimension `d ∈ 1:h.N`. Violations throw an
-`AssertionError`. Pass `check_args=false` to skip these checks in trusted
-inner loops.
+With `check_args=true` (the default), the method checks if the event `t` satisfies
+`h.tmin <= t < h.tmax`, if it occurs after the last existing event time for that
+dimention, and (if `h` is multivariate) if the dimension lie in `d ∈ 1:h.N`.
+Pass `check_args=false` to skip these checks in trusted inner loops.
 """
 function Base.push!(h::History, t::Real, m=nothing, d=nothing; check_args=true)
     if check_args
-        @assert h.tmin <= t < h.tmax
-        @assert (length(h) == 0) || (h.times[end] <= t)
-        @assert (d === nothing && h.N == 1) || 1 <= d <= h.N
+        if !(h.tmin <= t < h.tmax)
+            throw(
+                DomainError(
+                    (t, h.tmin, h.tmax),
+                    "Event time must lie in the half-open interval [tmin, tmax).",
+                ),
+            )
+        end
+        if !((d === nothing && h.N == 1) || (d isa Integer && 1 <= d <= h.N))
+            throw(
+                DomainError(
+                    d,
+                    "Event dimension must be between 1 and h.N for a multivariate history, or omitted for a univariate history.",
+                ),
+            )
+        end
+        if !(isempty(h) || event_times(h, d)[end] < t)
+            throw(
+                DomainError(
+                    (t, event_times(h, d)[end]),
+                    "Event time must be strictly greater than the previous event time in the same dimension.",
+                ),
+            )
+        end
     end
     push!(h.times, t)
     push!(h.marks, m)
diff --git a/src/mark_distributions.jl b/src/mark_distributions.jl
new file mode 100644
index 0000000..8174fcd
--- /dev/null
+++ b/src/mark_distributions.jl
@@ -0,0 +1,68 @@
+# Type definition
+"Abstract type for defining mark distributions not in `Distributions.jl`"
+abstract type AbstractMarkDistribution end
+
+const PointProcessMarkDistribution = Union{Distribution,AbstractMarkDistribution}
+
+## Standard implementations. Override as needed
+"""
+    mark_distribution(md, t, h)
+
+Compute the distribution of marks at time `t` after history `h`.
+
+Remark: This method must work for empty histories.
+"""
+function mark_distribution end
+
+function mark_distribution(md::AbstractMarkDistribution, t, h::History)
+    return error(
+        "Type $(typeof(md)) subtypes `AbstractMarkDistribution` but has " *
+        "not implemented the required `mark_distribution(md, t, h)` method.",
+    )
+end
+
+"""
+    sample_mark(rng, md, t, h)
+
+Return one sample from the distribution of marks at time `t` after history `h`, using the random number generator `rng`.
+"""
+function sample_mark(rng::AbstractRNG, md::PointProcessMarkDistribution, t, h::History)
+    return rand(rng, mark_distribution(md, t, h))
+end
+
+function sample_mark(md::PointProcessMarkDistribution, t, h::History)
+    return sample_mark(default_rng(), md, t, h)
+end
+
+"The type of the marks returned by the mark distribution"
+function Base.eltype(md::AbstractMarkDistribution)
+    return typeof(sample_mark(md, 0.0, History(0.0, 1.0)))
+end
+
+"The likelihood of a mark `m` occurring in an event at time `t` after history `h`"
+function DensityInterface.densityof(md::PointProcessMarkDistribution, t, h::History, m)
+    return densityof(mark_distribution(md, t, h), m)
+end
+
+# Support for `Distributions.jl`
+mark_distribution(d::Distribution, t, h::History) = d
+
+StatsAPI.fit(D::Type{<:Distribution}, h::History) = fit(D, h.marks)
+
+# Struct for non-marked processes
+"Mark distribution for non-marked processes. Always return the mark `nothing`."
+struct NoMarks <: AbstractMarkDistribution end
+
+mark_distribution(::NoMarks, t, h::History) = Dirac(nothing)
+
+sample_mark(::AbstractRNG, ::NoMarks, t, ::History) = nothing
+
+Base.eltype(::NoMarks) = Nothing
+
+StatsAPI.fit(::Type{NoMarks}, h) = NoMarks()
+
+StatsAPI.fit(::Type{NoMarks}, marks, weights) = NoMarks()
+
+DensityInterface.densityof(::NoMarks, t, ::History, ::Nothing) = 1.0
+
+DensityInterface.densityof(::NoMarks, t, ::History, m) = 0.0
diff --git a/src/multivariate/independent_multivariate.jl b/src/multivariate/independent_multivariate.jl
index b3bbdd4..95e4bc6 100644
--- a/src/multivariate/independent_multivariate.jl
+++ b/src/multivariate/independent_multivariate.jl
@@ -19,7 +19,7 @@ Base.ndims(pp::IndependentMultivariateProcess) = length(pp.processes)
 
 ## AbstractPointProcess interface
 function ground_intensity(pp::IndependentMultivariateProcess, t, h, d)
-    return ground_intensity(pp.processes[d], t, h)
+    return ground_intensity(pp.processes[d], t, History(h, d))
 end
 
 function ground_intensity(pp::IndependentMultivariateProcess, t, h)
@@ -27,7 +27,7 @@ function ground_intensity(pp::IndependentMultivariateProcess, t, h)
 end
 
 function mark_distribution(pp::IndependentMultivariateProcess, t, h, d)
-    return mark_distribution(pp.processes[d], t, h)
+    return mark_distribution(pp.processes[d], t, History(h, d))
 end
 
 function mark_distribution(pp::IndependentMultivariateProcess, t, h)
@@ -39,7 +39,7 @@ function mark_distribution(pp::IndependentMultivariateProcess, t)
 end
 
 function intensity(pp::IndependentMultivariateProcess, m, t, h, d)
-    return intensity(pp.processes[d], m, t, h)
+    return intensity(pp.processes[d], m, t, History(h, d))
 end
 
 function intensity(pp::IndependentMultivariateProcess, m, t, h)
@@ -49,10 +49,11 @@ end
 function log_intensity(pp::IndependentMultivariateProcess, m, t, h, d)
     return log(intensity(pp, m, t, h, d))
 end
+
 log_intensity(pp::IndependentMultivariateProcess, m, t, h) = log.(intensity(pp, m, t, h))
 
 function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h, d)
-    return ground_intensity_bound(pp.processes[d], t, h)
+    return ground_intensity_bound(pp.processes[d], t, History(h, d))
 end
 
 function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h)
@@ -60,7 +61,7 @@ function ground_intensity_bound(pp::IndependentMultivariateProcess, t, h)
 end
 
 function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b, d)
-    return integrated_ground_intensity(pp.processes[d], h, a, b)
+    return integrated_ground_intensity(pp.processes[d], History(h, d), a, b)
 end
 
 function integrated_ground_intensity(pp::IndependentMultivariateProcess, h, a, b)
@@ -76,10 +77,7 @@ function DensityInterface.logdensityof(pp::IndependentMultivariateProcess, h::Hi
 end
 
 function time_change(h::History{R,M}, pp::IndependentMultivariateProcess) where {R<:Real,M}
-    histories = [
-        time_change(History(event_times(h, d), h.tmin, h.tmax), pp.processes[d]) for
-        d in 1:ndims(pp)
-    ]
+    histories = [time_change(History(h, d), pp.processes[d]) for d in 1:ndims(pp)]
     tmax = maximum(histories[d].tmax for d in 1:ndims(pp))
     return History(
         [histories[d].times for d in 1:ndims(pp)],
@@ -101,11 +99,7 @@ end
 
 function StatsAPI.fit(proc_types::Vector, h::History{R,M}; kwargs...) where {R<:Real,M}
     processes = [
-        fit(
-            proc_types[d],
-            History(event_times(h, d), h.tmin, h.tmax, event_marks(h, d));
-            kwargs...,
-        ) for d in eachindex(proc_types)
+        fit(proc_types[d], History(h, d); kwargs...) for d in eachindex(proc_types)
     ]
     return IndependentMultivariateProcess(processes)
 end
diff --git a/src/multivariate/multivariate_poisson_process.jl b/src/multivariate/multivariate_poisson_process.jl
index 43a2190..052d3db 100644
--- a/src/multivariate/multivariate_poisson_process.jl
+++ b/src/multivariate/multivariate_poisson_process.jl
@@ -16,7 +16,7 @@ function PoissonProcess(
 end
 
 function PoissonProcess(λ::AbstractVector{R}; check_args::Bool=true) where {R<:Real}
-    return PoissonProcess(λ, [Dirac(nothing) for _ in eachindex(λ)]; check_args=check_args)
+    return PoissonProcess(λ, [NoMarks() for _ in eachindex(λ)]; check_args=check_args)
 end
 
 function PoissonProcess(
diff --git a/src/simulation.jl b/src/simulation.jl
index 08d4fd5..dfed8ca 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -15,27 +15,24 @@ end
     simulate_ogata(rng, pp, tmin, tmax)
 
 Simulate a temporal point process `pp` on interval `[tmin, tmax)` using Ogata's algorithm.
-
-# Technical Remark
-To infer the type of the marks, the implementation assumes that there is method of `mark_distribution` without the argument `h` such that it corresponds to the distribution of marks in case the history is empty.
 """
 function simulate_ogata(
-    rng::AbstractRNG, pp::AbstractPointProcess, tmin::T, tmax::T
+    rng::AbstractRNG, pp::AbstractUnivariateProcess, tmin::T, tmax::T
 ) where {T<:Real}
-    M = typeof(rand(mark_distribution(pp, tmin)))
+    M = eltype(pp.mark_dist)
     h = History(; times=T[], marks=M[], tmin=tmin, tmax=tmax)
     t = tmin
     while t < tmax
         B, L = ground_intensity_bound(pp, t + eps(t), h)
-        τ = B > 0 ? rand(rng, Exponential(inv(B))) : typemax(inv(B))
+        τ = B > 0 ? T(rand(rng, Exponential(inv(B)))) : typemax(inv(B))
         if τ > L
             t = t + L
         elseif τ <= L
             U_max = ground_intensity(pp, t + τ, h) / B
             U = rand(rng, typeof(U_max))
             if U < U_max
-                m = rand(rng, mark_distribution(pp, t + τ, h))
                 if t + τ < tmax
+                    m = sample_mark(pp.mark_dist, t + τ, h)
                     push!(h, t + τ, m; check_args=false)
                 end
             end
@@ -45,15 +42,6 @@ function simulate_ogata(
     return h
 end
 
-"""
-    simulate([rng,] pp, tmin, tmax)
-
-Alias for `simulate_ogata`.
-"""
-function simulate(rng::AbstractRNG, pp::AbstractPointProcess, tmin, tmax)
-    return simulate_ogata(rng, pp, tmin, tmax)
-end
-
 function simulate(pp::AbstractPointProcess, args...; kwargs...)
     return simulate(default_rng(), pp, args...; kwargs...)
 end
diff --git a/src/univariate/poisson/fit.jl b/src/univariate/poisson/fit.jl
index 0c34b0e..52ade53 100644
--- a/src/univariate/poisson/fit.jl
+++ b/src/univariate/poisson/fit.jl
@@ -1,16 +1,5 @@
 ## Fit MLE
 
-#=
-A separate `fit` method for unmarked Poisson processes is needed, because
-`Distributions.jl` does not provide a `fit` method for the `Dirac` distribution
-=#
-function StatsAPI.fit(
-    ::Type{PoissonProcess{R,Dirac{Nothing}}}, ss::PoissonProcessStats{R1,R2}; kwargs...
-) where {R<:Real,R1<:Real,R2<:Real}
-    λ = convert(R, ss.nb_events / ss.duration)
-    return PoissonProcess(λ, Dirac(nothing))
-end
-
 function StatsAPI.fit(
     ::Type{PoissonProcess{R,D}}, ss::PoissonProcessStats{R1,R2}; kwargs...
 ) where {R<:Real,D,R1<:Real,R2<:Real}
diff --git a/src/univariate/poisson/inhomogeneous/fit.jl b/src/univariate/poisson/inhomogeneous/fit.jl
index 9c4418c..024a652 100644
--- a/src/univariate/poisson/inhomogeneous/fit.jl
+++ b/src/univariate/poisson/inhomogeneous/fit.jl
@@ -82,18 +82,6 @@ function StatsAPI.fit(
     return from_params(F, optimal_params; intensity_kwargs...)
 end
 
-function StatsAPI.fit(
-    ::Type{<:InhomogeneousPoissonProcess{F,Dirac{Nothing}}},
-    h::History,
-    init_params;
-    integration_config=IntegrationConfig(),
-    kwargs...,
-) where {F<:ParametricIntensity}
-    # For unmarked processes, skip fitting the mark distribution
-    intensity = fit(F, h, init_params; integration_config=integration_config, kwargs...)
-    return InhomogeneousPoissonProcess(intensity, Dirac(nothing))
-end
-
 function StatsAPI.fit(
     ::Type{<:InhomogeneousPoissonProcess{F,D}},
     h::History,
@@ -102,7 +90,7 @@ function StatsAPI.fit(
     kwargs...,
 ) where {F<:ParametricIntensity,D}
     # Fit mark distribution independently
-    mark_dist = fit(D, h.marks)
+    mark_dist = fit(D, h)
 
     intensity = fit(F, h, init_params; integration_config=integration_config, kwargs...)
     return InhomogeneousPoissonProcess(intensity, mark_dist)
@@ -128,32 +116,6 @@ function StatsAPI.fit(
     return fit(pptype, combined, args...)
 end
 
-function StatsAPI.fit(
-    ::Type{<:InhomogeneousPoissonProcess{PiecewiseConstantIntensity{R},Dirac{Nothing}}},
-    h::History,
-    breakpoints::Vector;
-    kwargs...,
-) where {R}
-    # Count events in each bin
-    n_bins = length(breakpoints) - 1
-    rates = Vector{R}(undef, n_bins)
-    for i in 1:n_bins
-        # Count events in bin [breakpoints[i], breakpoints[i+1])
-        bin_start = breakpoints[i]
-        bin_end = breakpoints[i + 1]
-        bin_width = bin_end - bin_start
-
-        # Count events in this bin
-        count = sum(bin_start .<= h.times .< bin_end)
-
-        # Estimate rate as count / duration
-        rates[i] = count / bin_width
-    end
-
-    intensity = PiecewiseConstantIntensity(breakpoints, rates)
-    return InhomogeneousPoissonProcess(intensity, Dirac(nothing))
-end
-
 function StatsAPI.fit(
     ::Type{<:InhomogeneousPoissonProcess{PiecewiseConstantIntensity{R},D}},
     h::History,
@@ -161,7 +123,7 @@ function StatsAPI.fit(
     kwargs...,
 ) where {R,D}
     # Fit mark distribution independently
-    mark_dist = fit(D, h.marks)
+    mark_dist = fit(D, h)
 
     # Count events in each bin
     n_bins = length(breakpoints) - 1
diff --git a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
index 254bfa3..e922db4 100644
--- a/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
+++ b/src/univariate/poisson/inhomogeneous/inhomogeneous_poisson_process.jl
@@ -50,9 +50,7 @@ end
 function InhomogeneousPoissonProcess(
     f::F; integration_config::C=IntegrationConfig()
 ) where {F,C}
-    return InhomogeneousPoissonProcess{F,Dirac{Nothing},C}(
-        f, Dirac(nothing), integration_config
-    )
+    return InhomogeneousPoissonProcess{F,NoMarks,C}(f, NoMarks(), integration_config)
 end
 
 function Base.show(io::IO, pp::InhomogeneousPoissonProcess)
@@ -61,21 +59,12 @@ function Base.show(io::IO, pp::InhomogeneousPoissonProcess)
     )
 end
 
-## Access methods
-
-"""
-Mark distribution is time-independent for this implementation.
-"""
-mark_distribution(pp::InhomogeneousPoissonProcess) = pp.mark_dist
-
 ## AbstractPointProcess interface
 
 ground_intensity(pp::InhomogeneousPoissonProcess, t, h) = pp.intensity_function(t)
-mark_distribution(pp::InhomogeneousPoissonProcess, t, h) = pp.mark_dist
-mark_distribution(pp::InhomogeneousPoissonProcess, t) = pp.mark_dist
 
 function intensity(pp::InhomogeneousPoissonProcess, m, t, h)
-    return ground_intensity(pp, t, h) * densityof(mark_distribution(pp, t, h), m)
+    return ground_intensity(pp, t, h) * densityof(pp.mark_dist, t, h, m)
 end
 
 """
diff --git a/src/univariate/poisson/poisson_process.jl b/src/univariate/poisson/poisson_process.jl
index f6896cf..acd54a3 100644
--- a/src/univariate/poisson/poisson_process.jl
+++ b/src/univariate/poisson/poisson_process.jl
@@ -12,7 +12,7 @@ Homogeneous temporal Poisson process with arbitrary mark distribution.
 
     PoissonProcess(λ, mark_dist)
 """
-struct PoissonProcess{R<:Real,D} <: AbstractUnivariateProcess
+struct PoissonProcess{R<:Real,D<:PointProcessMarkDistribution} <: AbstractUnivariateProcess
     λ::R
     mark_dist::D
 
@@ -29,7 +29,7 @@ struct PoissonProcess{R<:Real,D} <: AbstractUnivariateProcess
 end
 
 function PoissonProcess(λ::R; check_args::Bool=true) where {R<:Real}
-    return PoissonProcess(λ, Dirac(nothing); check_args=check_args)
+    return PoissonProcess(λ, NoMarks(); check_args=check_args)
 end
 
 PoissonProcess() = PoissonProcess(1.0)
@@ -39,40 +39,23 @@ function Base.show(io::IO, pp::PoissonProcess)
 end
 
 ## Access
-ground_intensity(pp::PoissonProcess) = pp.λ
-mark_distribution(pp::PoissonProcess) = pp.mark_dist
-
-## Intensity functions
-function intensity(pp::PoissonProcess, m)
-    return ground_intensity(pp) * densityof(mark_distribution(pp), m)
-end
-
-function log_intensity(pp::PoissonProcess, m)
-    return log(ground_intensity(pp)) + logdensityof(mark_distribution(pp), m)
-end
+ground_intensity(pp::PoissonProcess, t, h) = pp.λ
 
 ## Time change
 function time_change(h::History, pp::PoissonProcess)
-    times = (h.times .- h.tmin) .* ground_intensity(pp)
-    tmax = (h.tmax - h.tmin) * ground_intensity(pp)
+    times = (h.times .- h.tmin) .* pp.λ
+    tmax = (h.tmax - h.tmin) * pp.λ
     return History(; times=times, tmin=0, tmax=tmax, marks=h.marks)
 end
 
 ## Implementing AbstractPointProcess interface
-
-ground_intensity(pp::PoissonProcess, t, h) = ground_intensity(pp)
-mark_distribution(pp::PoissonProcess, t, h) = mark_distribution(pp)
-mark_distribution(pp::PoissonProcess, t) = mark_distribution(pp) # For simulate_ogata
-intensity(pp::PoissonProcess, m, t, h) = intensity(pp, m)
-log_intensity(pp::PoissonProcess, m, t, h) = log_intensity(pp, m)
-
 function ground_intensity_bound(pp::PoissonProcess{R}, t::T, h) where {R,T<:Real}
     U = promote_type(R, T)
-    B = U(ground_intensity(pp))
+    B = U(pp.λ)
     L = typemax(U)
     return (B, L)
 end
 
 function integrated_ground_intensity(pp::PoissonProcess, h, a, b)
-    return ground_intensity(pp) * (b - a)
+    return pp.λ * (b - a)
 end
diff --git a/src/univariate/poisson/simulation.jl b/src/univariate/poisson/simulation.jl
index aa22d8d..bf4872e 100644
--- a/src/univariate/poisson/simulation.jl
+++ b/src/univariate/poisson/simulation.jl
@@ -1,5 +1,11 @@
 function simulate(rng::AbstractRNG, pp::PoissonProcess, tmin::T, tmax::T) where {T<:Real}
-    times = simulate_poisson_times(rng, ground_intensity(pp), tmin, tmax)
-    marks = [rand(rng, mark_distribution(pp)) for _ in 1:length(times)]
-    return History(; times=times, marks=marks, tmin=tmin, tmax=tmax)
+    h = History(T[], tmin, tmax, eltype(pp.mark_dist)[])
+    inter_dist = Exponential(inv(pp.λ))
+    t = T(rand(rng, inter_dist))
+    while t < tmax
+        m = sample_mark(pp.mark_dist, t, h)
+        push!(h, t, m; check_args=false)
+        t += rand(rng, inter_dist)
+    end
+    return h
 end
diff --git a/test/bounded_point_process.jl b/test/bounded_point_process.jl
index 494d756..718b9ae 100644
--- a/test/bounded_point_process.jl
+++ b/test/bounded_point_process.jl
@@ -8,9 +8,26 @@ h = simulate(rng, bpp)
 @test max_time(bpp) == 1000.0
 @test ground_intensity(bpp, 0, h) ≈ intensities
 @test mark_distribution(bpp, 100.0, h, 1) == Normal()
-@test mark_distribution(bpp, 0.0) == fill(Normal(), length(intensities))
+@test mark_distribution(bpp, 0.0, h) == fill(Normal(), length(intensities))
 @test intensity(bpp, 0.0, 0.0, h, 1) ≈ pdf(Normal(), 0.0) * intensities[1]
 @test log_intensity(bpp, 1.0, 1.0, h, 2) ≈ log(pdf(Normal(), 1.0)) + log(intensities[2])
 
 @test all(ground_intensity_bound(bpp, 243.0, h, 1) .≈ (intensities[1], Inf))
 @test integrated_ground_intensity(bpp, h, 342, 598) ≈ intensities * (598 - 342)
+
+struct FakePoisson <: AbstractUnivariateProcess
+    λ::Float64
+    mark_dist::NoMarks
+end
+
+PointProcesses.ground_intensity(fp::FakePoisson, t, h::History) = fp.λ
+function PointProcesses.integrated_ground_intensity(fp::FakePoisson, h::History, a, b)
+    return fp.λ * (b - a)
+end
+PointProcesses.ground_intensity_bound(fp::FakePoisson, t, h::History) = (fp.λ, Inf)
+
+pp = BoundedPointProcess(FakePoisson(1.0, NoMarks()), 0.0, 10.0)
+h = History([1.0, 2.0], 0.0, 10.0)
+
+@test simulate(pp) isa History
+@test time_change(h, pp).times == h.times
diff --git a/test/history.jl b/test/history.jl
index 5ed977d..339dc5e 100644
--- a/test/history.jl
+++ b/test/history.jl
@@ -1,6 +1,20 @@
 @testset "Univariate History" begin
+    # Constructors
+    h_empty1 = History(0.0, 1.0)
+    h_empty2 = History(0.0, 1.0, 2)
+
+    @test h_empty1 isa History{Float64,Any}
+    @test h_empty2 isa History{Float64,Any}
+    @test ndims(h_empty1) == 1
+    @test ndims(h_empty2) == 2
+    @test isempty(h_empty1)
+    @test isempty(h_empty2)
+
     h = History([0.2, 0.8, 1.1], 0.0, 2.0, ["a", "b", "c"])
 
+    @test h isa History{Float64,String}
+
+    # Access
     @test duration(h) == 2.0
     @test nb_events(h) == 3
     @test nb_events(h, 1.0, 2.0) == 1
@@ -11,20 +25,23 @@
     @test event_times(h) == [0.2, 0.8, 1.1]
     @test event_times(h, 0.2, 0.8) == [0.2]
     @test event_times(h, 0.8, 0.2) == []
+    @test event_times(h, nothing) == h.times
     @test event_marks(h) == ["a", "b", "c"]
     @test event_marks(h, 0.2, 0.8) == ["a"]
     @test event_marks(h, 0.8, 0.2) == []
     @test ndims(h) == 1
     @test event_dims(h) == fill(nothing, 3)
 
+    # Interface
+    @test_throws DomainError push!(h, 100.0, "d")
+    @test_throws DomainError push!(h, 100.0, "d", 2)
+    @test_throws DomainError push!(h, 1.0, "d", 2)
     push!(h, 1.7, "d")
 
     @test has_events(h, 1.5, 2.0)
 
     h2 = History(; times=[2.3], marks=["e"], tmin=2.0, tmax=2.5)
-    h3 = History(; times=[[1], [2, 2.5]], tmin=0, tmax=3)
 
-    @test (event_marks(h3) == fill(nothing, 3)) && (event_dims(h3) == [1, 2, 2])
     @test string(h2) == "History{Float64,String} with 1 events on interval [2.0, 2.5)"
 
     h_cat = cat(h, h2)
@@ -65,9 +82,24 @@ end
     @test event_marks(h_multi) == ["a", "c", "b", "d"]
     @test event_dims(h_multi) == [1, 2, 1, 2]
 
+    h_multi1 = History(h_multi, 1)
+    h_multi2 = History(h_multi, 2)
+
+    @test event_times(h_multi1) == times1
+    @test event_times(h_multi2) == times2
+    @test event_marks(h_multi1) == marks1
+    @test event_marks(h_multi2) == marks2
+
     @test_throws DomainError History(rand(3), 0, 1, rand(3), [1, 2, 3], 2)
     @test event_dims(History([[0.5]], 0, 1)) == [nothing]
 
+    @test_throws DomainError History(
+        [1.0, 1.0, 2.0, 3.0, 4.0], 0.0, 5.0, fill(nothing, 5), [1, 1, 1, 2, 1], 2
+    )
+    @test_throws DomainError History(
+        [1.0, 1.0, 1.0, 2.0, 3.0], 0.0, 5.0, fill(nothing, 5), [1, 2, 1, 2, 1], 2
+    )
+
     # Test dimension-specific methods
     @test event_times(h_multi, 1) == [0.1, 0.5]
     @test event_times(h_multi, 2) == [0.2, 0.8]
@@ -83,6 +115,7 @@ end
     @test event_dims(h_multi, 0.0, 0.3) == [1, 2]
 
     # Test push! with dimension
+    @test_throws DomainError push!(h_multi, 1.0, "d", nothing)
     push!(h_multi, 0.85, "e", 1)
     @test nb_events(h_multi, 1) == 3
     @test event_marks(h_multi, 1) == ["a", "b", "e"]
diff --git a/test/hypothesis_tests.jl b/test/hypothesis_tests.jl
index c980245..95b4367 100644
--- a/test/hypothesis_tests.jl
+++ b/test/hypothesis_tests.jl
@@ -1,17 +1,16 @@
-h1 = History([1, 2, 3, 4], 0, 5)
-h_empty = History(Float64[], 0, 2)
-PP = PoissonProcess{Float32,Dirac{Nothing}}
+h1 = History([1.0, 2.0, 3.0, 4.0], 0.0, 5.0)
+h_empty = History(Float64[], 0.0, 2.0)
+PP = PoissonProcess{Float32,NoMarks}
 pp = PoissonProcess()
 
 @testset "Statistics" begin
     @test statistic(KSDistance{Uniform}, pp, h1) ≈ 0.2
     @test statistic(KSDistance{Exponential}, pp, h1) ≈ 1 - exp(-1)
-
     @test statistic(KSDistance{Uniform}, pp, h_empty) ≈ 1
     @test statistic(KSDistance{Exponential}, pp, h_empty) ≈ 1
 end
 
-h2 = History(collect(0:999), 0, 1000)
+h2 = History(collect(0.0:999.0), 0.0, 1000.0)
 
 @testset "BootstrapTest" begin
     @test_throws ArgumentError BootstrapTest(KSDistance{Uniform}, PP, h_empty)
diff --git a/test/inhomogeneous_poisson_process.jl b/test/inhomogeneous_poisson_process.jl
index 70077d7..0f9fc80 100644
--- a/test/inhomogeneous_poisson_process.jl
+++ b/test/inhomogeneous_poisson_process.jl
@@ -334,7 +334,6 @@ end
         @test ground_intensity(pp, 0.0, h_empty) ≈ 1.0
         @test ground_intensity(pp, 2.0, h_empty) ≈ 2.0
         @test mark_distribution(pp, 0.0, h_empty) isa Normal
-        @test mark_distribution(pp) isa Normal
     end
 
     @testset "Simulation" begin
@@ -1110,8 +1109,7 @@ end
 
         @test pp isa InhomogeneousPoissonProcess
         @test pp.intensity_function === intensity
-        @test pp.mark_dist isa Dirac{Nothing}
-        @test pp.mark_dist.value === nothing
+        @test pp.mark_dist isa NoMarks
         @test pp.integration_config isa IntegrationConfig
 
         # Test that it simulates correctly
@@ -1127,7 +1125,7 @@ end
 
         @test pp isa InhomogeneousPoissonProcess
         @test pp.intensity_function === intensity
-        @test pp.mark_dist isa Dirac{Nothing}
+        @test pp.mark_dist isa NoMarks
         @test pp.integration_config === custom_config
         @test pp.integration_config.abstol == 1e-10
         @test pp.integration_config.reltol == 1e-10
@@ -1156,27 +1154,26 @@ end
     end
 end
 
-@testset "Fitting - PiecewiseConstant with Dirac (unmarked)" begin
+@testset "Fitting - PiecewiseConstant with no marks" begin
     @testset "Fit with explicit breakpoints vector" begin
         # Generate data from known piecewise constant process
         breakpoints_true = [0.0, 3.0, 7.0, 10.0]
         rates_true = [1.5, 4.0, 2.5]
         intensity_true = PiecewiseConstantIntensity(breakpoints_true, rates_true)
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 10.0)
 
         # Fit using same breakpoints
         pp_est = fit(
-            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks},
             h,
             breakpoints_true,
         )
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function isa PiecewiseConstantIntensity
-        @test pp_est.mark_dist isa Dirac{Nothing}
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
         @test length(pp_est.intensity_function.rates) == 3
 
         # Check breakpoints are preserved
@@ -1196,13 +1193,13 @@ end
         # Create uneven breakpoints
         custom_breakpoints = [0.0, 2.0, 5.0, 12.0]
         intensity_true = PiecewiseConstantIntensity([0.0, 6.0, 12.0], [3.0, 2.0])
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 12.0)
 
         # Fit with different breakpoints than data generation
         pp_est = fit(
-            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks},
             h,
             custom_breakpoints,
         )
@@ -1210,7 +1207,7 @@ end
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function.breakpoints == custom_breakpoints
         @test length(pp_est.intensity_function.rates) == 3
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # All rates should be positive
         @test all(r >= 0 for r in pp_est.intensity_function.rates)
@@ -1223,7 +1220,7 @@ end
         h = History([4.5, 4.8, 5.2, 5.5, 5.9], 0.0, 10.0)
 
         pp_est = fit(
-            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks},
             h,
             breakpoints,
         )
@@ -1243,19 +1240,19 @@ end
         # Edge case: single bin (constant rate)
         breakpoints = [0.0, 10.0]
         intensity_true = PolynomialIntensity([3.0])  # constant
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 10.0)
 
         pp_est = fit(
-            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float64},NoMarks},
             h,
             breakpoints,
         )
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test length(pp_est.intensity_function.rates) == 1
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # Rate should be approximately count / duration
         expected_rate = length(h.times) / 10.0
@@ -1268,7 +1265,7 @@ end
         h = History([1.0, 2.0, 6.0, 7.0, 8.0], 0.0, 10.0)
 
         pp_est = fit(
-            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float32},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PiecewiseConstantIntensity{Float32},NoMarks},
             h,
             breakpoints_f32,
         )
@@ -1278,25 +1275,22 @@ end
     end
 end
 
-@testset "Fitting - ParametricIntensity with Dirac (unmarked)" begin
-    @testset "Polynomial intensity with Dirac marks" begin
+@testset "Fitting - ParametricIntensity with no marks" begin
+    @testset "Polynomial intensity with no marks" begin
         # Generate data from linear process
         intensity_true = PolynomialIntensity([2.0, 0.3])
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 20.0)
 
         # Fit polynomial intensity
         pp_est = fit(
-            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}},
-            h,
-            [2.0, 0.3],
+            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks}, h, [2.0, 0.3]
         )
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function isa PolynomialIntensity
-        @test pp_est.mark_dist isa Dirac{Nothing}
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # Check parameters are reasonable
         @test abs(pp_est.intensity_function.coefficients[1] - 2.0) < 1.0
@@ -1311,13 +1305,13 @@ end
     @testset "Polynomial intensity with log link" begin
         # Generate data with log link (ensures positivity)
         intensity_true = PolynomialIntensity([1.5, 0.05]; link=:log)
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 30.0)
 
         # Fit with log link
         pp_est = fit(
-            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks},
             h,
             [1.5, 0.05];
             link=:log,
@@ -1325,7 +1319,7 @@ end
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function.link === :log
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # Verify positivity over range
         @test all(pp_est.intensity_function(t) > 0 for t in 0.0:0.5:30.0)
@@ -1335,24 +1329,23 @@ end
         @test abs(pp_est.intensity_function.coefficients[2] - 0.05) < 0.1
     end
 
-    @testset "Exponential intensity with Dirac marks" begin
+    @testset "Exponential intensity with no marks" begin
         # Generate data from exponential process
         intensity_true = ExponentialIntensity(3.0, 0.08)
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 15.0)
 
         # Fit exponential intensity
         pp_est = fit(
-            InhomogeneousPoissonProcess{ExponentialIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{ExponentialIntensity{Float64},NoMarks},
             h,
             [log(3.0), 0.08],
         )
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function isa ExponentialIntensity
-        @test pp_est.mark_dist isa Dirac{Nothing}
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # Check parameters are in reasonable range
         @test 0.5 * intensity_true.a <=
@@ -1364,16 +1357,16 @@ end
         @test all(pp_est.intensity_function(t) > 0 for t in 0.0:0.5:15.0)
     end
 
-    @testset "Sinusoidal intensity with Dirac marks" begin
+    @testset "Sinusoidal intensity with no marks" begin
         # Generate data from sinusoidal process
         intensity_true = SinusoidalIntensity(6.0, 2.5, 2π, 0.0)
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 8.0)
 
         # Fit sinusoidal intensity
         pp_est = fit(
-            InhomogeneousPoissonProcess{SinusoidalIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{SinusoidalIntensity{Float64},NoMarks},
             h,
             [log(6.0), 0.4, 0.0];
             ω=2π,
@@ -1381,8 +1374,7 @@ end
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function isa SinusoidalIntensity
-        @test pp_est.mark_dist isa Dirac{Nothing}
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
         @test pp_est.intensity_function.ω ≈ 2π
 
         # Check constraint a >= |b| is satisfied
@@ -1399,44 +1391,44 @@ end
     @testset "Custom integration config" begin
         # Test that integration_config is properly passed through
         intensity_true = PolynomialIntensity([2.5, 0.2])
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 20.0)
 
         # Fit with custom integration config
         custom_config = IntegrationConfig(; abstol=1e-10, reltol=1e-10, maxiters=5000)
         pp_est = fit(
-            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks},
             h,
             [2.5, 0.2];
             integration_config=custom_config,
         )
 
         @test pp_est isa InhomogeneousPoissonProcess
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # Verify fitting succeeded with custom config
         @test pp_est.intensity_function isa PolynomialIntensity
         @test abs(pp_est.intensity_function.coefficients[1] - 2.5) < 1.5
     end
 
-    @testset "Quadratic polynomial with Dirac marks" begin
+    @testset "Quadratic polynomial with no marks" begin
         # Test higher-degree polynomial
         intensity_true = PolynomialIntensity([1.0, 0.5, 0.05])
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 15.0)
 
         # Fit quadratic
         pp_est = fit(
-            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks},
             h,
             [1.0, 0.5, 0.05],
         )
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test length(pp_est.intensity_function.coefficients) == 3
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
 
         # Check all coefficients are in reasonable range
         for i in 1:3
@@ -1449,13 +1441,13 @@ end
     @testset "Pass through additional kwargs" begin
         # Test that optimizer kwargs are properly passed through
         intensity_true = PolynomialIntensity([2.0, 0.1])
-        pp_true = InhomogeneousPoissonProcess(intensity_true, Dirac(nothing))
+        pp_true = InhomogeneousPoissonProcess(intensity_true, NoMarks())
 
         h = simulate(rng, pp_true, 0.0, 25.0)
 
         # Fit with custom optimizer settings
         pp_est = fit(
-            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},Dirac{Nothing}},
+            InhomogeneousPoissonProcess{PolynomialIntensity{Float64},NoMarks},
             h,
             [2.0, 0.1];
             optimizer=LBFGS(),
@@ -1464,7 +1456,7 @@ end
 
         @test pp_est isa InhomogeneousPoissonProcess
         @test pp_est.intensity_function isa PolynomialIntensity
-        @test pp_est.mark_dist.value === nothing
+        @test pp_est.mark_dist isa NoMarks
     end
 end
 
diff --git a/test/mark_distributions.jl b/test/mark_distributions.jl
new file mode 100644
index 0000000..f0abd21
--- /dev/null
+++ b/test/mark_distributions.jl
@@ -0,0 +1,45 @@
+@testset "Distributions.jl" begin
+    h = History(sort(rand(10)), 0, 1, rand(10))
+    md = Normal()
+    @test md isa PointProcessMarkDistribution
+    @test mark_distribution(md, 0.0, h) == Normal()
+    @test fit(Normal, h) isa Normal
+    @test fit(Exponential, h) isa Exponential
+    @test densityof(md, 0.0, h, 0.0) == densityof(md, 0.0)
+    sample1 = sample_mark(Random.seed!(1), md, 0.0, h)
+    sample2 = rand(Random.seed!(1), Normal())
+    @test sample1 == sample2
+end
+
+@testset "NoMarks" begin
+    h = History(sort(rand(10)), 0, 1, rand(10))
+    md = NoMarks()
+    @test md isa AbstractMarkDistribution
+    @test md isa PointProcessMarkDistribution
+    @test mark_distribution(md, 0.0, h) isa Dirac{Nothing}
+    @test fit(NoMarks, h) == NoMarks()
+    @test fit(NoMarks, rand(10), rand(10)) == NoMarks()
+    @test densityof(md, 0.0, h, 0.0) == 0.0
+    @test densityof(md, 0.0, h, nothing) == 1.0
+    @test sample_mark(md, 0.0, h) === nothing
+    @test eltype(NoMarks()) == Nothing
+end
+
+@testset "Custom Mark Distribution" begin
+    struct TestMark <: AbstractMarkDistribution
+        t::Float64
+    end
+
+    h = History(0.0, 1.0)
+    md = TestMark(1.0)
+
+    @test_throws "not implemented" mark_distribution(md, 0.0, h)
+    @test_throws "not implemented" sample_mark(md, 0.0, h)
+    @test_throws "not implemented" eltype(md)
+    @test_throws "not implemented" densityof(md, 0.0, h, 0.0)
+
+    PointProcesses.mark_distribution(md::TestMark, t, h::History) = Normal(md.t)
+    @test mark_distribution(md, 0.0, h) == Normal(1.0)
+    @test eltype(md) == Float64
+    @test densityof(md, 0.0, h, 0.0) == densityof(Normal(1.0), 0.0)
+end
diff --git a/test/multivariate_poisson_process.jl b/test/multivariate_poisson_process.jl
index c0c6a58..416843d 100644
--- a/test/multivariate_poisson_process.jl
+++ b/test/multivariate_poisson_process.jl
@@ -37,17 +37,17 @@ end
     @test integrated_ground_intensity(pp1, h, 0, 1) == [1.0, 2.0, 3.0]
     @test integrated_ground_intensity(pp1, h, 0, 1, 1) == 1.0
 
-    h1 = simulate(rng, pp, 0.0, 1000.0)
+    h1 = simulate(rng, pp2, 0.0, 1000.0)
 
     f1(λ) = logdensityof(PoissonProcess(λ), h1)
-    gf = ForwardDiff.gradient(f1, 3 * ones(10))
+    gf = ForwardDiff.gradient(f1, 10 * ones(10))
 
     @test all(gf .< 0)
 end
 
 @testset "Simulation" begin
     pp0 = PoissonProcess([0.0, 1.0, 0.0])
-    bpp = BoundedPointProcess(pp1, 0, 1000)
+    bpp = BoundedPointProcess(pp1, 0.0, 1000.0)
     h1 = simulate(rng, pp1, 0.0, 1000.0)
     h2 = simulate(rng, pp0, 0.0, 1000.0)
     h3 = simulate(rng, bpp)
@@ -67,22 +67,29 @@ end
 end
 
 @testset "Fitting" begin
-    h1 = simulate(rng, pp1, 0.0, 1000.0)
-    h2 = History(fill(rand(100) .* 1000, 3), 0, 1000, fill(rand(100), 3))
+    h1 = simulate(rng, pp2, 0.0, 1000.0)
+    h2 = History(fill(rand(100) .* 1000, 3), 0.0, 1000.0)
 
-    pp_est1 = fit(fill(PoissonProcess{Float64,Normal}, 3), h1)
-    pp_est2 = fit(fill(PoissonProcess{Float64,Normal}, 3), h2)
+    pp_est1 = fit(fill(PoissonProcess{Float64,NoMarks}, 3), h1)
+    pp_est2 = fit(fill(PoissonProcess{Float64,NoMarks}, 3), h2)
     λ_est1 = [pp_est1.processes[d].λ for d in 1:ndims(pp_est1)]
     λ_est2 = [pp_est2.processes[d].λ for d in 1:ndims(pp_est2)]
 
     λ_error1 = mean(abs, λ_est1 - λ)
     λ_error2 = mean(abs, λ_est2 - λ)
 
-    l = logdensityof(pp, h1)
+    l = logdensityof(pp2, h1)
     l_est = logdensityof(pp_est1, h1)
 
     @test λ_error1 < λ_error2
     @test l_est > l
+
+    pp_marks = PoissonProcess([1, 1], Normal())
+    h_marks = simulate(pp_marks, 0.0, 1000.0)
+    pp_est_marks = fit(fill(PoissonProcess{Float64,Normal}, 2), h_marks)
+    @test pp_est_marks isa MultivariatePoissonProcess
+    @test pp_est_marks.processes[1].mark_dist isa Normal
+    @test logdensityof(pp_est_marks, h_marks) >= logdensityof(pp_marks, h_marks)
 end
 
 @testset "Time change" begin
diff --git a/test/poisson_process.jl b/test/poisson_process.jl
index 53fcdb7..4b994b6 100644
--- a/test/poisson_process.jl
+++ b/test/poisson_process.jl
@@ -13,8 +13,8 @@ rng = Random.seed!(63)
     @test pp2.λ == 1.0
     @test pp3.λ == 1.0
     @test pp1.mark_dist == Normal()
-    @test pp2.mark_dist == Dirac(nothing)
-    @test pp3.mark_dist == Dirac(nothing)
+    @test pp2.mark_dist == NoMarks()
+    @test pp3.mark_dist == NoMarks()
     @test_throws DomainError PoissonProcess(-1.0)
 end
 
@@ -27,7 +27,7 @@ pp = PoissonProcess(1.0, Normal())
     h = History([0.1, 0.5], 0, 1, [0.0, 0.0])
 
     @test ground_intensity(pp, 0, h) == 1.0
-    @test mark_distribution(pp) == Normal()
+    @test mark_distribution(pp, 0.0, h) == Normal()
     @test intensity(pp, 0.5, 0, h) == 1.0 * pdf(Normal(), 0.5)
     @test log_intensity(pp, 0.5, 0, h) == log(1.0) + logpdf(Normal(), 0.5)
     @test integrated_ground_intensity(pp, h, 0.0, 1.0) == 1.0
diff --git a/test/runtests.jl b/test/runtests.jl
index 1a47516..4c0d52e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -44,6 +44,9 @@ DocMeta.setdocmeta!(PointProcesses, :DocTestSetup, :(using PointProcesses); recu
     @testset verbose = true "History" begin
         include("history.jl")
     end
+    @testset verbose = false "MarkDistributions" begin
+        include("mark_distributions.jl")
+    end
     @testset verbose = true "Bounded" begin
         include("bounded_point_process.jl")
     end