Commit 8c9e8b30 authored by Javier's avatar Javier

Merge branch 'alejandro'

parents 7c0e8824 46a49854
......@@ -11,7 +11,7 @@ include("juobs_obs.jl")
export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!
export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs, hess_reduce, uwcholesky, transpose, tridiag_reduction, make_positive_def, invert_covar_matrix
export corr_obs, corr_obs_TSM, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine
export corr_obs, corr_obs_TSM, corr_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine, bayesian_av
export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
end # module
......@@ -473,6 +473,306 @@ 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, pl::Bool, data::Bool; 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, pl::Bool, data::Bool; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
bayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Array{Int64}, pl::Bool, data::Bool; 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 of them a weight. The function saves the `first` fit parameter of your function, and then it does the weighted average of it 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. Then, for each fit interval choice, the function explores the two fit functions. This means that for each fit interval choice you get two results: one for the first fit funcction, and another for the second. You can also use a vector of functions and a vector of k (numer of parameters of each funtion) to apply the bayesian averaging method to multiple functions.
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. If `data` is `true`, then returns 4 objects: weighted average, systematic error, a vector with the results of the fit for each fit interval choice, and a vector with the weights associated to each fit.
```@example
@.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, data, weights) = bayesian_av(fun,x,tmin_array,tmax_array,k,pl=true,data=true)
@.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)
@.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, pl::Bool=false, data::Bool=false; 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}()
mods = Array{String,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 = gen_chisq(fun,x,dy)
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)
isnothing(wpm) ? uwerr(up[1]) : 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])
push!(mods,string("[", INDEX+1, ",", j, "]"))
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)
if pl
x = 1:length(p1)
y = value.(p1)
dy = err.(p1)
v = value(p1_mean)
e = err(p1_mean)
figure()
fill_between(1:length(p1), v-e, v+e, color="green", alpha=0.75)
errorbar(mods, y, dy, fmt="x", color="black")
ylabel(L"$p_1$")
xlabel(L"model")
display(gcf())
figure()
errorbar(mods, weight_model, 0*dy, color="green")
ylabel(L"$weight$")
xlabel(L"model")
display(gcf())
end
if !data
return (p1_mean, systematic_err)
else
return (p1_mean, systematic_err, p1, weight_model)
end
end
function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k1::Int64, k2::Int64, pl::Bool=false, data::Bool=false; 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}()
mods = Array{String,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 = gen_chisq(fun1,x,dy)
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)
isnothing(wpm) ? uwerr(up[1]) : 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])
push!(mods,string("[", INDEX+1, ",", j, "]"))
p00 = [0.5 for i in 1:1:k2]
chisq = gen_chisq(fun2,x,dy)
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])
push!(mods,string("[", INDEX+1, ",", j, "]"))
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)
if pl
x = 1:length(p1)
y = value.(p1)
dy = err.(p1)
v = value(p1_mean)
e = err(p1_mean)
figure()
fill_between(1:length(p1), v-e, v+e, color="green", alpha=0.75)
errorbar(mods, y, dy, fmt="x", color="black")
ylabel(L"$p_1$")
xlabel(L"model")
display(gcf())
figure()
errorbar(mods, weight_model, 0*dy, color="green")
ylabel(L"$weight$")
xlabel(L"model")
display(gcf())
end
if !data
return (p1_mean, systematic_err)
else
return (p1_mean, systematic_err, p1, weight_model)
end
end
function bayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Array{Int64}, pl::Bool=false, data::Bool=false; 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}()
mods = Array{String,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
for indice in 1:length(fun)
p00 = [0.5 for i in 1:1:k[indice]]
chisq = gen_chisq(fun[indice],x,dy)
fit = curve_fit(fun[indice],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)
isnothing(wpm) ? uwerr(up[1]) : uwerr(up[1],wpm)
chi2 = sum(fit.resid.^2) * dof(fit) / chi_exp
push!(AIC, chi2 + 2*k[indice] + 2*Ncut)
push!(chi2chi2exp, chi2 / dof(fit))
push!(p1, up[1])
push!(mods,string("[", INDEX+1, ",", j, "]"))
end
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)
if pl
x = 1:length(p1)
y = value.(p1)
dy = err.(p1)
v = value(p1_mean)
e = err(p1_mean)
figure()
fill_between(1:length(p1), v-e, v+e, color="green", alpha=0.75)
errorbar(mods, y, dy, fmt="x", color="black")
ylabel(L"$p_1$")
xlabel(L"model")
display(gcf())
figure()
errorbar(mods, weight_model, 0*dy, color="green")
ylabel(L"$weight$")
xlabel(L"model")
display(gcf())
end
if !data
return (p1_mean, systematic_err)
else
return (p1_mean, systematic_err, p1, weight_model)
end
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