Commit d0004abc by AlejandroSaezGonzalvo

### Update juobs_tools.jl

`Bayesian averaging method supports now averaging over two fit functions`
parent cfa5a31f
 ... ... @@ -309,19 +309,30 @@ 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) bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::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. You can also use as input two fit functions, and two values of `k`, one for each function. 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 @.fun(x,p) = p[1] * x ^0 k = 1 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) @.fun1(x,p) = p[1] * x ^0 @.fun2(x,p) = p[1] + p[2] * exp( - p[3] * (x)) k1 = 1 k2 = 3 tmin_array = [10,11,12,13,14,15] tmax_array = [80,81,82,83,84,85] (average, systematics) = bayesian_av(fun1,fun2,x,tmin_array,tmax_array,k1,k2) ``` """ 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) ... ... @@ -374,6 +385,68 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, end function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::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:k1] chisq = fit_defs(fun1,x,W) fit = curve_fit(fun1,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*k1 + 2*Ncut) push!(chi2chi2exp, chi2 / dof(fit)) push!(p1, up[1]) p00 = [0.5 for i in 1:1:k2] chisq = fit_defs(fun2,x,W) fit = curve_fit(fun2,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*k2 + 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) #; uwerr(systematic_err) 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) ... ...
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment