Update in bayesian_av to allow plots

src/juobs_tools.jl

@@ -307,16 +307,22 @@ function plat_av(obs::Vector{uwreal}, plat::Vector{Int64}, wpm::Union{Dict{Int64
 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(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::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; 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; 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.
+
+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.
+
 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
@@ -324,7 +330,7 @@ The method returns two objects: first, the weighted average as an uwreal object,
 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)
+(average, systematics) = bayesian_av(fun,x,tmin_array,tmax_array,k,pl=true)
 
 @.fun1(x,p) = p[1] * x ^0
 @.fun2(x,p) = p[1] + p[2] * exp( - p[3] * (x))
@@ -335,7 +341,7 @@ 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)
+function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64}, tmax_array::Array{Int64}, k::Int64, pl::Bool=false; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 
     weight_model = Array{Float64,1}()
     AIC = Array{Float64,1}()
@@ -381,11 +387,27 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
     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(x, y, dy, fmt="x", color="black")
+        ylabel(L"$p_1$")
+        xlabel(L"model")
+    
+        display(gcf())
+    end
+
     return (p1_mean, systematic_err)  
 
 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)
+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; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
 
     weight_model = Array{Float64,1}()
     AIC = Array{Float64,1}()
@@ -440,7 +462,23 @@ function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_arra
     
     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)
+    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(x, y, dy, fmt="x", color="black")
+        ylabel(L"$p_1$")
+        xlabel(L"model")
+    
+        display(gcf())
+    end
 
     return (p1_mean, systematic_err)