Update juobs_tools.jl

Bayesian average method for one fit function added to juobs_tools.jl

src/juobs_tools.jl

@@ -306,6 +306,74 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64
     return av 
 end
 
+@doc raw"""
+    bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
+Computes bayesian average of data. For a given fit function, it explores choices of fit intervals, assigning each a weight. Then it does the weighted average and assigns a systematic. See https://arxiv.org/abs/2008.01069 
+The function takes as input the fit intervals to explore. 
+`tmin_array` is an array of integers with the lower bounds on the fit intervals to explore, ***ordered from lower to higher***.
+`tmax_array` is an array of integers with the upper bounds on the fit intervals to explore, ***ordered from lower to higher***.
+`k` is the number of parameters of the fit function to use.
+The method returns two objects: first, the weighted average as an uwreal object, with mean value and statistichal error. The second object returned is the systematic error coming from the fit interval variation.
+
+```@example
+@.fun(x,p) = p[1] + p[2] * exp( - p[3] * (x)) + p[4] * exp( - p[5] * (dim - x)) 
+k = 5
+tmin_array = [10,11,12,13,14,15]
+tmax_array = [80,81,82,83,84,85]
+(average, systematics) = bayesian_av(fun,x,tmin_array,tmax_array,k)
+```
+"""
+function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+
+    weight_model = Array{Float64,1}()
+    AIC = Array{Float64,1}()
+    chi2chi2exp = Array{Float64,1}()
+    p1 = Array{uwreal,1}()
+
+    if tmax_array[end] > length(y)
+        error("Error: upper bound for the fits is bigger than last data point")
+    end
+
+    total = length(y)
+    isnothing(wpm) ? uwerr.(y) : for i in 1:length(y) uwerr(y[i],wpm) end
+    
+    for INDEX in tmin_array ## vary tmin
+        for j in tmax_array ## vary tmax
+            try
+                x = [i for i in INDEX+1:1:j]	
+	            yy = y[INDEX+1:1:j]
+	            Ncut = total - length(x)
+                dy = err.(yy)
+                W = 1 ./ dy .^2
+
+                p00 = [0.5 for i in 1:1:k]
+                chisq = fit_defs(fun,x,W)
+                fit = curve_fit(fun,x,value.(yy),W,p00)
+                isnothing(wpm) ? (up,chi_exp) = fit_error(chisq,coef(fit),yy) : (up,chi_exp) = fit_error(chisq,coef(fit),yy,wpm)
+                uwerr(up[1],wpm)
+                chi2 = sum(fit.resid.^2) * dof(fit) / chi_exp
+                push!(AIC, chi2 + 2*k + 2*Ncut)
+                push!(chi2chi2exp, chi2 / dof(fit))
+                push!(p1, up[1])
+            catch e
+                @warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point")
+            end
+        end
+    end
+
+    offset = minimum(AIC)
+    AIC = AIC .- offset
+    weight_model = exp.(-0.5 .* AIC)
+    
+    p1_mean = sum(p1 .* weight_model)/sum(weight_model) ; isnothing(wpm) ? uwerr(p1_mean) : uwerr(p1_mean,wpm) 
+    weight_model = weight_model ./ sum(weight_model)
+    systematic_err = sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2) 
+
+    return (p1_mean, systematic_err)  
+
+end
+
 function lin_fit(x::Vector{<:Real}, v::Vector{Float64}, e::Vector{Float64})
     sig2 = e .* e
     S = sum(1 ./ sig2)