Skip to content

Commit

Permalink
Add explicit proposal_sigma based on both aveuncert and data spacing
Browse files Browse the repository at this point in the history
  • Loading branch information
@brenhinkeller
brenhinkeller committed May 27, 2024
1 parent 1b76661 commit 44748e7
Showing 1 changed file with 42 additions and 19 deletions.
61 changes: 42 additions & 19 deletions src/StratMetropolis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,18 @@
# Stratigraphic age constraints
Age = copy(smpl.Age)::Vector{Float64}
Age_sigma = copy(smpl.Age_sigma)::Vector{Float64}
aveuncert = nanmean(Age_sigma)
Height = copy(smpl.Height)::Vector{Float64}
Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
Chronometer = smpl.Chronometer
(bottom, top) = extrema(Height)
model_heights = bottom:resolution:top

aveuncert = nanmean(Age_sigma)
absdiff = diff(sort!(Age[Age_Sidedness.==0]))
maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)

if bounding>0
# If bounding is requested, add extrapolated top and bottom bounds to avoid
# issues with the stratigraphic markov chain wandering off to +/- infinity
Expand All @@ -66,7 +70,7 @@

# Run the Markov chain
ages = Normal.(Age, Age_sigma)
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, aveuncert, burnin, nsteps, sieve, Chronometer, systematic)
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)

# Crop the result
agedist = agedist[active_height_t,:]
Expand All @@ -89,13 +93,17 @@
# Stratigraphic age constraints. Type assertions for stability
Age = copy(smpl.Age)::Vector{Float64}
Age_sigma = copy(smpl.Age_sigma)::Vector{Float64}
aveuncert = nanmean(Age_sigma)
Height = copy(smpl.Height)::Vector{Float64}
Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
(bottom, top) = extrema(Height)
model_heights = bottom:resolution:top

aveuncert = nanmean(Age_sigma)
absdiff = diff(sort!(Age[Age_Sidedness.==0]))
maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)

if bounding>0
# If bounding is requested, add extrapolated top and bottom bounds to avoid
# issues with the stratigraphic markov chain wandering off to +/- infinity
Expand All @@ -118,7 +126,7 @@

# Run the Markov chain
ages = Normal.(Age, Age_sigma)
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, aveuncert, burnin, nsteps, sieve)
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)

# Crop the result
agedist = agedist[active_height_t,:]
Expand Down Expand Up @@ -165,7 +173,6 @@
# Stratigraphic age constraints
Age = copy(smpl.Age)::Vector{Float64}
Age_sigma = copy(smpl.Age_sigma)::Vector{Float64}
aveuncert = nanmean(Age_sigma)
Height = copy(smpl.Height)::Vector{Float64}
Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
Expand All @@ -174,6 +181,10 @@
(bottom, top) = extrema(Height)
model_heights = bottom:resolution:top

aveuncert = nanmean(Age_sigma)
absdiff = diff(sort!(Age[Age_Sidedness.==0]))
maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)

if bounding>0
# If bounding is requested, add extrapolated top and bottom bounds to avoid
Expand Down Expand Up @@ -201,7 +212,7 @@

# Run the Markov chain
ages = BilinearExponential.(eachcol(p))
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, aveuncert, burnin, nsteps, sieve, Chronometer, systematic)
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)

# Crop the result
agedist = agedist[active_height_t,:]
Expand All @@ -224,14 +235,18 @@
# Stratigraphic age constraints
Age = copy(smpl.Age)::Vector{Float64}
Age_sigma = copy(smpl.Age_sigma)::Vector{Float64}
aveuncert = nanmean(Age_sigma)
Height = copy(smpl.Height)::Vector{Float64}
Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
p = copy(smpl.Params)::Matrix{Float64}
(bottom, top) = extrema(Height)
model_heights = bottom:resolution:top

aveuncert = nanmean(Age_sigma)
absdiff = diff(sort!(Age[Age_Sidedness.==0]))
maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)

if bounding>0
# If bounding is requested, add extrapolated top and bottom bounds to avoid
# issues with the stratigraphic markov chain wandering off to +/- infinity
Expand All @@ -257,7 +272,7 @@

# Run the Markov chain
ages = BilinearExponential.(eachcol(p))
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, aveuncert, burnin, nsteps, sieve)
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)

# Crop the result
agedist = agedist[active_height_t,:]
Expand Down Expand Up @@ -304,14 +319,18 @@
# Stratigraphic age constraints
Age = copy(smpl.Age)::Vector{Float64}
Age_sigma = copy(smpl.Age_sigma)::Vector{Float64}
aveuncert = nanmean(Age_sigma)
Height = copy(smpl.Height)::Vector{Float64}
Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
p = copy(smpl.Params)::Matrix{Float64}
(bottom, top) = extrema(Height)
model_heights = bottom:resolution:top

aveuncert = nanmean(Age_sigma)
absdiff = diff(sort!(Age[Age_Sidedness.==0]))
maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)

if bounding>0
# If bounding is requested, add extrapolated top and bottom bounds to avoid
# issues with the stratigraphic markov chain wandering off to +/- infinity
Expand All @@ -338,7 +357,7 @@

# Run the Markov chain
ages = Radiocarbon.(Age, Age_sigma, (collect(c) for c in eachcol(p)))
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, aveuncert, burnin, nsteps, sieve)
agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)

# Crop the result
agedist = agedist[active_height_t,:]
Expand All @@ -364,14 +383,18 @@
# Stratigraphic age constraints
Age = copy(smpl.Age)::Vector{Float64}
Age_sigma = copy(smpl.Age_sigma)::Vector{Float64}
aveuncert = nanmean(Age_sigma)
Height = copy(smpl.Height)::Vector{Float64}
Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
p = copy(smpl.Params)::Matrix{Float64}
(bottom, top) = extrema(Height)
model_heights = bottom:resolution:top

aveuncert = nanmean(Age_sigma)
absdiff = diff(sort!(Age[Age_Sidedness.==0]))
maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)

if bounding>0
# If bounding is requested, add extrapolated top and bottom bounds to avoid
# issues with the stratigraphic markov chain wandering off to +/- infinity
Expand All @@ -398,7 +421,7 @@

# Run the Markov chain
ages = Radiocarbon.(Age, Age_sigma, (collect(c) for c in eachcol(p)))
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, aveuncert, burnin, nsteps, sieve)
agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, Age_Sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)

# Crop the result
agedist = agedist[active_height_t,:]
Expand Down Expand Up @@ -443,7 +466,7 @@
return ages
end

function stratmetropolis(Height, Height_sigma, model_heights::AbstractRange, Age_Sidedness, ages, model_ages, aveuncert, burnin::Integer, nsteps::Integer, sieve::Integer, Chronometer=nothing, systematic=nothing)
function stratmetropolis(Height, Height_sigma, model_heights::AbstractRange, Age_Sidedness, ages, model_ages, proposal_sigma, burnin::Integer, nsteps::Integer, sieve::Integer, Chronometer=nothing, systematic=nothing)
resolution = step(model_heights)
npoints = length(model_heights)

Expand Down Expand Up @@ -494,7 +517,7 @@
end
else
# Adjust one point at a time then resolve conflicts
r = randn() * aveuncert # Generate a random adjustment
r = randn() * proposal_sigma # Generate a random adjustment
chosen_point = ceil(Int, rand() * npoints) # Pick a point
model_agesₚ[chosen_point] += r
#Resolve conflicts
Expand Down Expand Up @@ -569,7 +592,7 @@
end
else
# Adjust one point at a time then resolve conflicts
r = randn() * aveuncert # Generate a random adjustment
r = randn() * proposal_sigma # Generate a random adjustment
chosen_point = ceil(Int, rand() * npoints) # Pick a point
model_agesₚ[chosen_point] += r
#Resolve conflicts
Expand Down Expand Up @@ -623,7 +646,7 @@
return agedist, lldist
end

function stratmetropolis(hiatus::HiatusData, Height, Height_sigma, model_heights::AbstractRange, Age_Sidedness, ages, model_ages, aveuncert, burnin::Integer, nsteps::Integer, sieve::Integer, Chronometer=nothing, systematic=nothing)
function stratmetropolis(hiatus::HiatusData, Height, Height_sigma, model_heights::AbstractRange, Age_Sidedness, ages, model_ages, proposal_sigma, burnin::Integer, nsteps::Integer, sieve::Integer, Chronometer=nothing, systematic=nothing)
resolution = step(model_heights)
npoints = length(model_heights)

Expand Down Expand Up @@ -693,7 +716,7 @@
end
else
# Adjust one point at a time then resolve conflicts
r = randn() * aveuncert # Generate a random adjustment
r = randn() * proposal_sigma # Generate a random adjustment
chosen_point = ceil(Int, rand() * npoints) # Pick a point
model_agesₚ[chosen_point] += r
#Resolve conflicts
Expand Down Expand Up @@ -787,7 +810,7 @@
end
else
# Adjust one point at a time then resolve conflicts
r = randn() * aveuncert # Generate a random adjustment
r = randn() * proposal_sigma # Generate a random adjustment
chosen_point = ceil(Int, rand() * npoints) # Pick a point
model_agesₚ[chosen_point] += r
#Resolve conflicts
Expand Down Expand Up @@ -858,4 +881,4 @@
end
update!(pgrs,nsteps*sieve)
return agedist, lldist, hiatusdist
end
end

0 comments on commit 44748e7

Please sign in to comment.