Skip to content

Commit

Permalink
Fix hiatusdist output
Browse files Browse the repository at this point in the history
  • Loading branch information
brenhinkeller committed May 27, 2024
1 parent 51f2f10 commit b1457b0
Show file tree
Hide file tree
Showing 2 changed files with 17 additions and 22 deletions.
26 changes: 17 additions & 9 deletions src/StratMetropolis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@
model_heights = model_heights[active_height_t]
mdl = StratAgeModel(model_heights, agedist)

return mdl, agedist, lldist, hiatusdist
return mdl, agedist, hiatusdist, lldist
end

## --- Stratigraphic MCMC model without hiatus, with distribution LL # # # # # #
Expand Down Expand Up @@ -264,7 +264,7 @@
model_heights = model_heights[active_height_t]
mdl = StratAgeModel(model_heights, agedist)

return mdl, agedist, lldist, hiatusdist
return mdl, agedist, hiatusdist, lldist
end

## --- Stratigraphic MCMC model without hiatus, for radiocarbon ages # # # # # #
Expand Down Expand Up @@ -405,7 +405,7 @@
model_heights = model_heights[active_height_t]
mdl = StratAgeModel(model_heights, agedist)

return mdl, agedist, lldist, hiatusdist
return mdl, agedist, hiatusdist, lldist
end

## --- # Internals of the Markov chain
Expand Down Expand Up @@ -716,9 +716,13 @@
if model_agesₚ[h-1] == model_agesₚ[h]
n = findclosestunequal(model_agesₚ, h)
if n < h
model_agesₚ[n:h-1] .= model_agesₚ[n]
@inbounds for i = n:h-1
model_agesₚ[i] = model_agesₚ[n]
end
elseif n > h
model_agesₚ[h:n] .= model_agesₚ[n]
@inbounds for i = h:n
model_agesₚ[i] = model_agesₚ[n]
end
end
end
end
Expand All @@ -740,7 +744,7 @@

# Add log likelihood for hiatus duration
@. durationₚ = model_agesₚ[closest_hiatus_unique - 1] - model_agesₚ[closest_hiatus_unique]
llₚ += normcdf_ll!(Hiatus_duration, Hiatus_duration_sigma, durationₚ)
llₚ += normcdf_ll(Hiatus_duration, Hiatus_duration_sigma, durationₚ)

# Accept or reject proposal based on likelihood
if log(rand(Float64)) < (llₚ - ll)
Expand Down Expand Up @@ -809,9 +813,13 @@
if model_agesₚ[h-1] == model_agesₚ[h]
n = findclosestunequal(model_agesₚ, h)
if n < h
model_agesₚ[n:h-1] .= model_agesₚ[n]
@inbounds for i = n:h-1
model_agesₚ[i] = model_agesₚ[n]
end
elseif n > h
model_agesₚ[h:n] .= model_agesₚ[n]
@inbounds for i = h:n
model_agesₚ[i] = model_agesₚ[n]
end
end
end
end
Expand All @@ -832,7 +840,7 @@

# Add log likelihood for hiatus duration
@. durationₚ = model_agesₚ[closest_hiatus_unique - 1] - model_agesₚ[closest_hiatus_unique]
llₚ += normcdf_ll!(Hiatus_duration, Hiatus_duration_sigma, durationₚ)
llₚ += normcdf_ll(Hiatus_duration, Hiatus_duration_sigma, durationₚ)

# Accept or reject proposal based on likelihood
if log(rand(Float64)) < (llₚ - ll)
Expand Down
13 changes: 0 additions & 13 deletions src/foo.jl

This file was deleted.

0 comments on commit b1457b0

Please sign in to comment.