Skip to content

Commit 8b185b3

Browse files
brenhinkellerbrenhinkeller
authored
Merge pull request #30 from brenhinkeller/cdf
Add methods for slower but more accurate cdf-based handling of one-sided age constraints
2 parents bd16159 + 46992e5 commit 8b185b3

File tree

6 files changed

+209
-148
lines changed

6 files changed

+209
-148
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "Chron"
22
uuid = "68885b1f-77b5-52a7-b2e7-6a8014c56b98"
33
authors = ["C. Brenhin Keller <cbkeller@dartmouth.edu>"]
4-
version = "0.5.5"
4+
version = "0.6.0"
55

66
[deps]
77
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

‎src/Objects.jl

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
inputSigmaLevel::Int
2222
Age_Unit::String
2323
Height_Unit::String
24+
Sidedness_Method::Symbol
2425
end
2526

2627
function ChronAgeData(nSamples::Integer)
@@ -43,6 +44,7 @@
4344
2, # i.e., are the data files 1-sigma or 2-sigma
4445
"Ma",
4546
"m",
47+
:cdf,
4648
)
4749
end
4850

@@ -55,6 +57,7 @@
5557
Chronometer::NTuple{N, Symbol}
5658
Age_Unit::String
5759
Height_Unit::String
60+
Sidedness_Method::Symbol
5861
end
5962

6063
function GeneralAgeData(nSamples::Integer)
@@ -67,9 +70,27 @@
6770
ntuple(i->:Chronometer, nSamples), # Age Types (e.g., :UPb or :ArAr)
6871
"Ma",
6972
"m",
73+
:cdf,
7074
)
7175
end
7276

77+
# A set of types to allow dispatch to different age sidedness methods
78+
abstract type Sidedness end
79+
struct CDFSidedness{T<:Real} <: Sidedness
80+
s::Vector{T}
81+
end
82+
struct FastSidedness{T<:Real} <: Sidedness
83+
s::Vector{T}
84+
end
85+
Base.getindex(A::Sidedness, args...) = getindex(A.s, args...)
86+
Base.setindex!(A::Sidedness, args...) = setindex!(A.s, args...)
87+
Base.eachindex(A::Sidedness) = eachindex(A.s)
88+
Base.size(A::Sidedness, args...) = size(A.s, args...)
89+
Base.length(A::Sidedness) = length(A.s)
90+
Base.ndims(A::Sidedness) = ndims(A.s)
91+
Base.ndims(::Type{<:Sidedness}) = 1
92+
Base.copyto!(A::Sidedness, args...) = copyto!(A.s, args...)
93+
7394
# One-sigma systematic uncertainty
7495
mutable struct SystematicUncertainty
7596
UPb::Float64

‎src/StratMetropolis.jl

Lines changed: 75 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -68,9 +68,16 @@
6868
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
6969
model_ages = a .+ b .* collect(model_heights)
7070

71+
# Select sidedness method
72+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
73+
FastSidedness(Age_Sidedness)
74+
else
75+
CDFSidedness(Age_Sidedness)
76+
end
77+
7178
# Run the Markov chain
7279
ages = Normal.(Age, Age_sigma)
73-
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
80+
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
7481

7582
# Crop the result
7683
agedist = agedist[active_height_t,:]
@@ -124,9 +131,16 @@
124131
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
125132
model_ages = a .+ b .* collect(model_heights)
126133

134+
# Select sidedness method
135+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
136+
FastSidedness(Age_Sidedness)
137+
else
138+
CDFSidedness(Age_Sidedness)
139+
end
140+
127141
# Run the Markov chain
128142
ages = Normal.(Age, Age_sigma)
129-
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
143+
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
130144

131145
# Crop the result
132146
agedist = agedist[active_height_t,:]
@@ -210,9 +224,16 @@
210224
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
211225
model_ages = a .+ b .* collect(model_heights)
212226

227+
# Select sidedness method
228+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
229+
FastSidedness(Age_Sidedness)
230+
else
231+
CDFSidedness(Age_Sidedness)
232+
end
233+
213234
# Run the Markov chain
214235
ages = BilinearExponential.(eachcol(p))
215-
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
236+
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
216237

217238
# Crop the result
218239
agedist = agedist[active_height_t,:]
@@ -270,9 +291,16 @@
270291
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
271292
model_ages = a .+ b .* collect(model_heights)
272293

294+
# Select sidedness method
295+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
296+
FastSidedness(Age_Sidedness)
297+
else
298+
CDFSidedness(Age_Sidedness)
299+
end
300+
273301
# Run the Markov chain
274302
ages = BilinearExponential.(eachcol(p))
275-
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
303+
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
276304

277305
# Crop the result
278306
agedist = agedist[active_height_t,:]
@@ -355,9 +383,16 @@
355383
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
356384
model_ages = a .+ b .* collect(model_heights)
357385

386+
# Select sidedness method
387+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
388+
FastSidedness(Age_Sidedness)
389+
else
390+
CDFSidedness(Age_Sidedness)
391+
end
392+
358393
# Run the Markov chain
359394
ages = Radiocarbon.(Age, Age_sigma, (collect(c) for c in eachcol(p)))
360-
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
395+
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
361396

362397
# Crop the result
363398
agedist = agedist[active_height_t,:]
@@ -419,9 +454,16 @@
419454
(a,b) = hcat(fill!(similar(Height), 1), Height) \ Age
420455
model_ages = a .+ b .* collect(model_heights)
421456

457+
# Select sidedness method
458+
sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
459+
FastSidedness(Age_Sidedness)
460+
else
461+
CDFSidedness(Age_Sidedness)
462+
end
463+
422464
# Run the Markov chain
423465
ages = Radiocarbon.(Age, Age_sigma, (collect(c) for c in eachcol(p)))
424-
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
466+
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
425467

426468
# Crop the result
427469
agedist = agedist[active_height_t,:]
@@ -438,17 +480,11 @@
438480
npoints = length(model_heights)
439481

440482
# Calculate log likelihood of initial proposal
441-
# Proposals younger than age constraint are given a pass if Age_Sidedness is -1 (maximum age)
442-
# proposals older than age constraint are given a pass if Age_Sidedness is +1 (minimum age)
483+
443484
sample_height = copy(Height)
444485
closest = findclosest(sample_height, model_heights)
445486
closest_model_ages = model_ages[closest]
446-
@inbounds for i ∈ eachindex(ages,closest_model_ages)
447-
if Age_Sidedness[i] == sign(closest_model_ages[i] - ages[i].μ)
448-
closest_model_ages[i] = ages[i].μ
449-
end
450-
end
451-
ll = strat_ll(closest_model_ages, ages)
487+
ll = strat_ll(closest_model_ages, ages, Age_Sidedness)
452488
ll += normpdf_ll(Height, Height_sigma, sample_height)
453489

454490
# Preallocate variables for MCMC proposals
@@ -473,14 +509,10 @@
473509

474510
if rand() < 0.1
475511
# Adjust heights
476-
@inbounds for i ∈ eachindex(sample_heightₚ)
512+
@inbounds for i ∈ eachindex(sample_heightₚ, closestₚ)
477513
sample_heightₚ[i] += randn() * Height_sigma[i]
478-
closestₚ[i] = round(Int,(sample_heightₚ[i] - model_heights[1])/resolution)+1
479-
if closestₚ[i] < 1 # Check we're still within bounds
480-
closestₚ[i] = 1
481-
elseif closestₚ[i] > npoints
482-
closestₚ[i] = npoints
483-
end
514+
closestₚ[i] = round(Int,(sample_heightₚ[i] - first(model_heights))/resolution)+1
515+
closestₚ[i] = max(min(closestₚ[i], lastindex(model_agesₚ)), firstindex(model_agesₚ))
484516
end
485517
else
486518
# Adjust one point at a time then resolve conflicts
@@ -505,16 +537,11 @@
505537

506538

507539
# Calculate log likelihood of proposal
508-
# Proposals younger than age constraint are given a pass if Age_Sidedness is -1 (maximum age)
509-
# proposal older than age constraint are given a pass if Age_Sidedness is +1 (minimum age)
510-
@inbounds for i ∈ eachindex(ages, closest_model_agesₚ)
540+
adjust!(agesₚ, Chronometer, systematic)
541+
@inbounds for i ∈ eachindex(closest_model_agesₚ, closestₚ)
511542
closest_model_agesₚ[i] = model_agesₚ[closestₚ[i]]
512-
if Age_Sidedness[i] == sign(closest_model_agesₚ[i] - ages[i].μ)
513-
closest_model_agesₚ[i] = ages[i].μ
514-
end
515543
end
516-
adjust!(agesₚ, Chronometer, systematic)
517-
llₚ = strat_ll(closest_model_agesₚ, agesₚ)
544+
llₚ = strat_ll(closest_model_agesₚ, agesₚ, Age_Sidedness)
518545
llₚ += normpdf_ll(Height, Height_sigma, sample_heightₚ)
519546

520547
# Accept or reject proposal based on likelihood
@@ -549,14 +576,10 @@
549576

550577
if rand() < 0.1
551578
# Adjust heights
552-
@inbounds for i ∈ eachindex(sample_heightₚ)
579+
@inbounds for i ∈ eachindex(sample_heightₚ, closestₚ)
553580
sample_heightₚ[i] += randn() * Height_sigma[i]
554-
closestₚ[i] = round(Int,(sample_heightₚ[i] - model_heights[1])/resolution)+1
555-
if closestₚ[i] < 1 # Check we're still within bounds
556-
closestₚ[i] = 1
557-
elseif closestₚ[i] > npoints
558-
closestₚ[i] = npoints
559-
end
581+
closestₚ[i] = round(Int,(sample_heightₚ[i] - first(model_heights))/resolution)+1
582+
closestₚ[i] = max(min(closestₚ[i], lastindex(model_agesₚ)), firstindex(model_agesₚ))
560583
end
561584
else
562585
# Adjust one point at a time then resolve conflicts
@@ -580,16 +603,11 @@
580603
end
581604

582605
# Calculate log likelihood of proposal
583-
# Proposals younger than age constraint are given a pass if Age_Sidedness is -1 (maximum age)
584-
# proposal older than age constraint are given a pass if Age_Sidedness is +1 (minimum age)
585-
@inbounds for i ∈ eachindex(ages, closest_model_agesₚ)
606+
adjust!(agesₚ, Chronometer, systematic)
607+
@inbounds for i ∈ eachindex(closest_model_agesₚ, closestₚ)
586608
closest_model_agesₚ[i] = model_agesₚ[closestₚ[i]]
587-
if Age_Sidedness[i] == sign(closest_model_agesₚ[i] - ages[i].μ)
588-
closest_model_agesₚ[i] = ages[i].μ
589-
end
590609
end
591-
adjust!(agesₚ, Chronometer, systematic)
592-
llₚ = strat_ll(closest_model_agesₚ, agesₚ)
610+
llₚ = strat_ll(closest_model_agesₚ, agesₚ, Age_Sidedness)
593611
llₚ += normpdf_ll(Height, Height_sigma, sample_heightₚ)
594612

595613
# Accept or reject proposal based on likelihood
@@ -619,17 +637,10 @@
619637
npoints = length(model_heights)
620638

621639
# Calculate log likelihood of initial proposal
622-
# Proposals younger than age constraint are given a pass if Age_Sidedness is -1 (maximum age)
623-
# proposals older than age constraint are given a pass if Age_Sidedness is +1 (minimum age)
624640
sample_height = copy(Height)
625641
closest = findclosest(sample_height, model_heights)
626642
closest_model_ages = model_ages[closest]
627-
@inbounds for i ∈ eachindex(ages)
628-
if Age_Sidedness[i] == sign(closest_model_ages[i] - ages[i].μ)
629-
closest_model_ages[i] = ages[i].μ
630-
end
631-
end
632-
ll = strat_ll(closest_model_ages, ages)
643+
ll = strat_ll(closest_model_ages, ages, Age_Sidedness)
633644
ll += normpdf_ll(Height, Height_sigma, sample_height)
634645

635646
# Ensure there is only one effective hiatus at most for each height node
@@ -673,14 +684,10 @@
673684

674685
if rand() < 0.1
675686
# Adjust heights
676-
@inbounds for i ∈ eachindex(sample_heightₚ)
687+
@inbounds for i ∈ eachindex(sample_heightₚ, closestₚ)
677688
sample_heightₚ[i] += randn() * Height_sigma[i]
678-
closestₚ[i] = round(Int,(sample_heightₚ[i] - model_heights[1])/resolution)+1
679-
if closestₚ[i] < 1 # Check we're still within bounds
680-
closestₚ[i] = 1
681-
elseif closestₚ[i] > npoints
682-
closestₚ[i] = npoints
683-
end
689+
closestₚ[i] = round(Int,(sample_heightₚ[i] - first(model_heights))/resolution)+1
690+
closestₚ[i] = max(min(closestₚ[i], lastindex(model_agesₚ)), firstindex(model_agesₚ))
684691
end
685692
else
686693
# Adjust one point at a time then resolve conflicts
@@ -718,16 +725,11 @@
718725

719726

720727
# Calculate log likelihood of proposal
721-
# Proposals younger than age constraint are given a pass if Age_Sidedness is -1 (maximum age)
722-
# proposal older than age constraint are given a pass if Age_Sidedness is +1 (minimum age)
723-
@inbounds for i ∈ eachindex(ages)
728+
adjust!(agesₚ, Chronometer, systematic)
729+
@inbounds for i ∈ eachindex(closest_model_agesₚ, closestₚ)
724730
closest_model_agesₚ[i] = model_agesₚ[closestₚ[i]]
725-
if Age_Sidedness[i] == sign(closest_model_agesₚ[i] - ages[i].μ)
726-
closest_model_agesₚ[i] = ages[i].μ
727-
end
728731
end
729-
adjust!(agesₚ, Chronometer, systematic)
730-
llₚ = strat_ll(closest_model_agesₚ, agesₚ)
732+
llₚ = strat_ll(closest_model_agesₚ, agesₚ, Age_Sidedness)
731733
llₚ += normpdf_ll(Height, Height_sigma, sample_heightₚ)
732734

733735
# Add log likelihood for hiatus duration
@@ -767,14 +769,10 @@
767769

768770
if rand() < 0.1
769771
# Adjust heights
770-
@inbounds for i ∈ eachindex(sample_heightₚ)
772+
@inbounds for i ∈ eachindex(sample_heightₚ, closestₚ)
771773
sample_heightₚ[i] += randn() * Height_sigma[i]
772-
closestₚ[i] = round(Int,(sample_heightₚ[i] - model_heights[1])/resolution)+1
773-
if closestₚ[i] < 1 # Check we're still within bounds
774-
closestₚ[i] = 1
775-
elseif closestₚ[i] > npoints
776-
closestₚ[i] = npoints
777-
end
774+
closestₚ[i] = round(Int,(sample_heightₚ[i] - first(model_heights))/resolution)+1
775+
closestₚ[i] = max(min(closestₚ[i], lastindex(model_agesₚ)), firstindex(model_agesₚ))
778776
end
779777
else
780778
# Adjust one point at a time then resolve conflicts
@@ -811,16 +809,11 @@
811809
end
812810

813811
# Calculate log likelihood of proposal
814-
# Proposals younger than age constraint are given a pass if Age_Sidedness is -1 (maximum age)
815-
# proposal older than age constraint are given a pass if Age_Sidedness is +1 (minimum age)
816-
@inbounds for i ∈ eachindex(ages)
812+
adjust!(agesₚ, Chronometer, systematic)
813+
@inbounds for i ∈ eachindex(closest_model_agesₚ, closestₚ)
817814
closest_model_agesₚ[i] = model_agesₚ[closestₚ[i]]
818-
if Age_Sidedness[i] == sign(closest_model_agesₚ[i] - ages[i].μ)
819-
closest_model_agesₚ[i] = ages[i].μ
820-
end
821815
end
822-
adjust!(agesₚ, Chronometer, systematic)
823-
llₚ = strat_ll(closest_model_agesₚ, agesₚ)
816+
llₚ = strat_ll(closest_model_agesₚ, agesₚ, Age_Sidedness)
824817
llₚ += normpdf_ll(Height, Height_sigma, sample_heightₚ)
825818

826819
# Add log likelihood for hiatus duration

0 commit comments

Comments
 (0)