Fitting

The fit function can be used to fit a model to a recorded series with either the Whittle or debiased Whittle likelihoods using the interior point Newton method as implemented in Optim.jl.

First load the package:

using OceanWaveSpectralFitting

As a basic example, consider a univariate Gaussian process with JONSWAP spectral density function. We can simulate such a process with the following:

n = 2304
Δ = 1/1.28
nreps = 1
α = 0.7
ωₚ = 1.1
γ = 3.3
r = 5.0
ts = simulate_gp(JONSWAP{1}(α,ωₚ,γ,r),n,Δ,nreps)[1]

WhittleLikelihoodInference.TimeSeries{Float64, 1}([0.33119740817539817, -0.12853479070253782, -0.47246934279789865, -0.5301457875276546, 0.05187566181562843, 0.8416111982127756, 0.6486789646392398, -0.20731813527834136, -0.591153499962582, -0.5820249337051248  …  0.05974814830183698, 0.23437032265464003, 0.2770150046120959, 0.12166956544817194, 0.1206692066364476, -0.10908959486960154, -0.337704279986251, -0.27178406279100636, -0.06519866914116158, 0.3749690671942457], 0.78125)

In this case, we have simulated a Gaussian process with JONSWAP spectral density function with parameters α = 0.7, ωₚ = 1.1, γ = 3.3, r = 5.0 of length 2304 sampled every 1/1.28 seconds. This corresponds to 30 minutes of data. The function simulate_gp will simulate a vector of nreps series, which is why we recover the first of these to get one series.

We can now fit a model with the following:

x₀ = [α,ωₚ,γ,r] .+ 0.1
res = fit(ts,model=JONSWAP{1},x₀=x₀)

 * Status: success

 * Candidate solution
    Final objective value:     -1.168053e+04

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 2.10e-05 ≰ 1.0e-08

 * Work counters
    Seconds run:   3  (vs limit Inf)
    Iterations:    34
    f(x) calls:    140
    ∇f(x) calls:   140

Here x₀ is the vector of initial parameter guesses. The estimated parameter vector can be recovered by doing:

x̂ = res.minimizer

4-element Vector{Float64}:
 0.9844319608459207
 1.1178179709575218
 2.063265684874429
 5.358748141543505

The full example is:

using OceanWaveSpectralFitting
n = 2304
Δ = 1/1.28
nreps = 1
α = 0.7
ωₚ = 1.1
γ = 3.3
r = 5.0
ts = simulate_gp(JONSWAP{1}(α,ωₚ,γ,r),n,Δ,nreps)[1]
x₀ = [α,ωₚ,γ,r] .+ 0.1
res = fit(ts,model=JONSWAP{1},x₀=x₀)
x̂ = res.minimizer

4-element Vector{Float64}:
 0.9844319608459207
 1.1178179709575218
 2.063265684874429
 5.358748141543505

For more options, see below:

fit(ts,Δ;model::Type{<:TimeSeriesModel},x₀,lowerΩcutoff,upperΩcutoff,x_lowerbounds,x_upperbounds,method,taper)
fit(timeseries::TimeSeries;model::Type{<:TimeSeriesModel},x₀,lowerΩcutoff,upperΩcutoff,x_lowerbounds,x_upperbounds,method,taper)