Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change ARStep to reverse AR parameters #441

Closed
wants to merge 5 commits into from
Closed

Change ARStep to reverse AR parameters #441

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
4e5fe42
reverse damp parameters
SamuelBrand1 Aug 28, 2024
2b1411f
Add a theoretical test for AR(2)
SamuelBrand1 Aug 28, 2024
dfd9f6e
Add rev_damp_prior field to AR
SamuelBrand1 Aug 28, 2024
ba8a4b3
Update AR.jl
SamuelBrand1 Aug 28, 2024
4d6513b
fix plot with change in sampled quanitity from AR
SamuelBrand1 Aug 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions EpiAware/src/EpiLatentModels/models/AR.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,23 +81,26 @@ Generate a latent AR series.
damp_AR ~ latent_model.damp_prior
ϵ_t ~ filldist(Normal(), n - p)

ar = accumulate_scan(ARStep(damp_AR), ar_init, σ_AR * ϵ_t)
ar_step = ARStep(reverse(damp_AR))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For efficieny I don't think we want to do this at run time?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To avoid doing reverse? That would suggest reversing at the constructor level I think.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

e.g. we add a rev_damp_prior field to AR and have damp_AR ~ latent_model.rev_damp_prior. Obvs we could just reverse damp_prior but that might get confusing for users?

ar = accumulate_scan(ar_step, ar_init, σ_AR * ϵ_t)

return ar
end

@doc raw"
The autoregressive (AR) step function struct
The autoregressive (AR) step function struct.

Note that the AR parameters are stored in reverse order.
"
struct ARStep{D <: AbstractVector{<:Real}} <: AbstractAccumulationStep
damp_AR::D
rev_damp_AR::D
end

@doc raw"
The autoregressive (AR) step function for use with `accumulate_scan`.
"
function (ar::ARStep)(state, ϵ)
new_val = dot(ar.damp_AR, state) + ϵ
new_val = dot(ar.rev_damp_AR, state) + ϵ
new_state = vcat(state[2:end], new_val)
return new_state
end
31 changes: 30 additions & 1 deletion EpiAware/test/EpiLatentModels/models/AR.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ end
end
end

@testitem "Testing AR process against theoretical properties" begin
@testitem "Testing AR(1) process against theoretical properties" begin
using DynamicPPL, Turing
using HypothesisTests: ExactOneSampleKSTest, pvalue
using Distributions
Expand Down Expand Up @@ -88,3 +88,32 @@ end
samples, Normal(theoretical_mean, sqrt(theoretical_var))) |> pvalue
@test ks_test_pval > 1e-6
end

@testitem "Testing AR(2) process against theoretical properties" begin
using DynamicPPL, Turing
using HypothesisTests: ExactOneSampleKSTest, pvalue
using Distributions

ar_model = AR(Normal(), HalfNormal(0.1), Normal(), p = 2)
n = 10_000
damp = [0.8, 0.1]
σ_AR = 1.0

theoretical_mean = 0.0
theoretical_var = σ_AR^2 / (1 - damp[1]^2 - damp[2]^2 -
2 * (damp[1]^2 * damp[2] / (1 - damp[2])))
theoretical_1step_cov = theoretical_var * damp[1] / (1 - damp[2])
init_Σ = [theoretical_var theoretical_1step_cov; theoretical_1step_cov theoretical_var]

# Draw initial values from the stationary distribution so process starts at stationarity
ar_init = rand(MvNormal(zeros(2), init_Σ))

model = generate_latent(ar_model, n)
fixed_model = fix(model, (σ_AR = σ_AR, damp_AR = damp, ar_init = ar_init))
# Draw samples from the model
X = fixed_model()

ks_test_pval = ExactOneSampleKSTest(
X, Normal(theoretical_mean, sqrt(theoretical_var))) |> pvalue
@test ks_test_pval > 1e-6
end
Loading