Skip to content

Instantly share code, notes, and snippets.

@briochemc
Last active February 5, 2021 05:11
Show Gist options
  • Save briochemc/06a8629d88566bbaea5692c012011079 to your computer and use it in GitHub Desktop.
Save briochemc/06a8629d88566bbaea5692c012011079 to your computer and use it in GitHub Desktop.
### A Pluto.jl notebook ###
# v0.12.20
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
el
end
end
# ╔═╡ 0a180f60-66de-11eb-2a4a-7fb399216f09
begin
using Pkg
Pkg.activate(mktempdir())
end
# ╔═╡ 355ec972-66de-11eb-039b-fdc7fb07f42f
begin
Pkg.add("Plots")
Pkg.add("MonteCarloMeasurements")
Pkg.add("PlutoUI")
Pkg.add("SpecialFunctions")
Pkg.add("Distributions")
using Distributions
using Plots
using MonteCarloMeasurements
using PlutoUI
using SpecialFunctions
end
# ╔═╡ 531148ee-66de-11eb-2b95-95b74d044ecd
begin
@enum FeType dissolved anthropogenic mineral
δFe_range = Dict(
dissolved => (-0.65, -0.30),
anthropogenic => (-1.8 , -1.6),
mineral => (0.1 , 0.7)
)
md"""
# A more quantitative evaluation of the anthropogenic fraction
In the Pinedo et al (2020) paper, the ranges for δFe values are
- dissolved: $(δFe_range[dissolved])
- anthropogenic: $(δFe_range[anthropogenic])
- mineral: $(δFe_range[mineral])
and are sort of shown below (the y axis does not matter for now).
"""
end
# ╔═╡ 7b51b040-66df-11eb-3fcd-4343186ac585
let
plt = plot()
for t in instances(FeType)
vspan!(plt, collect(δFe_range[t]), label=t, c=Int(t), α=0.5)
end
xlabel!(plt, "δFe (‰)")
plt
end
# ╔═╡ ddf2aa6c-66e2-11eb-1a90-e901044ab5d8
begin
md"""
Naturally, the idea is that the dissolved (green) values are obtained by mixing the 2 end-members, anthropogenic (blue) and mineral (red) sources.
The equation in the Pinedo et al (2020) manuscript for the anthropogenic fraction is in the Methods section:
$$\delta\mathsf{Fe}_\mathsf{dissolved} = (1 - f) \, \delta\mathsf{Fe}_\mathsf{mineral} + f \, \delta\mathsf{Fe}_\mathsf{anthropogenic}$$
Note the notation changes: **The anthropogenic fraction** is now just $f$.
(And because fractions must sum to 1, I also set $f_\mathsf{mineral} = 1 - f$.)
So the equation for $f$ is
$$f = \frac{\delta\mathsf{Fe}_\mathsf{dissolved} - \delta\mathsf{Fe}_\mathsf{mineral}}{\delta\mathsf{Fe}_\mathsf{anthropogenic}- \delta\mathsf{Fe}_\mathsf{mineral}}$$
So let's use Julia and make some random variables of δFe to generate a distribution for $f$!
"""
end
# ╔═╡ e6c2ca56-66e1-11eb-2998-9931de836769
begin
md"""
But instead of using extremas, let's assume we have a Gaussian distribution, which mostly lies with the range. I wrote a little function for that, so you just have to pick a probability of the random variable being within the range, from 50 to 99.5%:
Probability of being in the range: $(@bind probs Slider(50:0.1:99.9, show_value=true))
"""
end
# ╔═╡ e11f3e18-6745-11eb-1d50-ad84d8390aae
begin
MC(t::FeType, prob::Number) = MC(δFe_range[t], prob)
#MC(x, prob::Number) = μ(x) ± σ(x, prob)
MC(x, prob::Number) = Particles(1000000, Uniform(x...))
IC(x, p) = quantile(x, [(1-p)/2, 1-(1-p)/2])
μ(x) = (x[1] + x[2]) / 2
σ(x, prob::Number) = (x[2] - x[1]) / (2sqrt(2) * erfinv(prob))
function f(probs)
MCδFe = Dict((t => MC(t, probs) for t in instances(FeType))...)
100(MCδFe[dissolved] - MCδFe[mineral]) / (MCδFe[anthropogenic] - MCδFe[mineral])
end
end
# ╔═╡ d143c98c-66e1-11eb-3ccb-4bc539223694
let
plt = plot()
for t in instances(FeType)
x = δFe_range[t]
vspan!(plt, collect(x), label=t, c=Int(t), α=0.5)
μt = round(mean(MC(t,probs)), sigdigits=2)
σt = round(std(MC(t,probs)), sigdigits=2)
annotate!(plt, μ(x), Int(t) + 7, text("$μt ± $σt"))
histogram!(plt, MC(x, probs/100), c=Int(t), α=0.5, label="", normalize=true, bins=-2.5:0.01:2, xticks=-2.5:0.5:2, ylims=(0,10))
end
xlabel!(plt, "δFe (‰)")
plt
end
# ╔═╡ ae1831a8-66e3-11eb-326b-f7a6cfbb5c3e
let
_f = f(probs/100)
μf, σf = mean(_f), std(_f)
bins = 0:1:100
plt = plot()
[(ic = IC(_f.particles, cf/100); vspan!(plt, ic, lab="$(cf)% confidence: $(round.(ic, sigdigits=3))", α=1.1-cf/100)) for cf in (68.27, 95.45, 99.73, 100.0)]
histogram!(plt, _f; title="Fraction of anthropogenic Fe is ($(round(μf, sigdigits=2)) ± $(round(σf, sigdigits=2)))%", xlab="%", label="probability", normalize=true, bins, xlims=extrema(bins), xticks=0:10:100, c=:gray, lc=:black)
plt
end
# ╔═╡ Cell order:
# ╟─0a180f60-66de-11eb-2a4a-7fb399216f09
# ╟─355ec972-66de-11eb-039b-fdc7fb07f42f
# ╟─531148ee-66de-11eb-2b95-95b74d044ecd
# ╟─7b51b040-66df-11eb-3fcd-4343186ac585
# ╟─ddf2aa6c-66e2-11eb-1a90-e901044ab5d8
# ╟─e6c2ca56-66e1-11eb-2998-9931de836769
# ╟─d143c98c-66e1-11eb-3ccb-4bc539223694
# ╟─ae1831a8-66e3-11eb-326b-f7a6cfbb5c3e
# ╟─e11f3e18-6745-11eb-1d50-ad84d8390aae
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment