Pinedo_etal_PNAS_2020_anthropogenic_fraction.jl

### 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