Skip to content

J
juobs
Commit f3ec800f authored by Alessandro 's avatar Alessandro
Browse files
Options

bayesian av implemented in my branch, new Corr constructor, minor modifications

parent d6bd07c8
Hide whitespace changes
Inline Side-by-side
Showing with 325 additions and 73 deletions
constants/juobs_const.jl
View file @ f3ec800f
......@@ -9,7 +9,7 @@ const ZM_data = [1.9684, 1.9935, 2.0253, 2.0630, 2.0814]
const ZM_error = [35, 27, 33, 38, 45] .* 1e-4
const ZM_tm_data = [2.6047, 2.6181, 2.6312, 2.6339, 2.6127]
const ZM_tm_error = [42, 33, 42, 48, 55] .* 1e-4
# 1808.09236
#1808.09236
const ZA_data = [0.75642, 0.76169, 0.76979, 0.78378, 0.79667]
const ZA_err = [72, 93, 43, 47, 47] .*1e-5
const ZV_data = [0.72259, 0.72845, 0.73886, 0.75574, 0.77074]
......@@ -18,16 +18,21 @@ const ZV_err = [74, 89, 46, 48, 49] .* 1e-5
const ZS_over_ZP_data = [1.3497, 1.2914, 1.2317, 1.1709, 1.1343]
const ZS_over_ZP_err = [83, 64, 48, 23, 25 ] .* 1e-4
#1608.08900
const b_values2 = [3.40, 3.46, 3.55, 3.70]
const t0_data = [2.86, 3.659, 5.164, 8.595]
const t0_error = [11, 16, 18, 29] .* 1e-3
# const b_values2 = [3.40, 3.46, 3.55, 3.70]
const t0_data = [2.86, 3.659, 5.164, 8.595, 14.040]
const t0_error = [11, 16, 18, 29, 49] .* 1e-3
const t0_ph_value = [0.413]
const t0_ph_error = ones(1,1) .* 5e-3
#2101.10969
const RM_data = [2.335, 1.869, 1.523, 1.267, 1.149]
const RM_err = [31, 19, 14, 16, 18] .* 1e-3
#1906.03445 LCP1 for heavy quarks
const BA_MINUS_BP_data = [-0.324, -0.265, -0.196, -0.119, -0.073]
const BA_MINUS_BP_err = [17, 14, 14, 14, 12] .*1e-3
const ZDATA = [0.8981, 0.9468, 1.0015, 1.0612, 1.0971]
const ZERR = [35, 35, 30, 17, 18 ] .* 1e-4
## taking into account correlations in zm_tm
#1802.05243 taking into account correlations in zm_tm
C = [[0.375369e-6, 0.429197e-6, -0.186896e-5] [0.429197e-6, 0.268393e-4, 0.686776e-4] [-0.186896e-5, 0.686776e-4, 0.212386e-3]]
z = cobs([0.348629, 0.020921, 0.070613], C, "z")
ZP(b) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
......@@ -43,12 +48,15 @@ zm_tm_covar(b) = Mrat / ZP(b)
#Cov matrices
const C1 = zeros(5, 5)
const C2 = zeros(5, 5)
const C3 = zeros(4, 4)
const C3 = zeros(5, 5)
const C4 = zeros(7, 7)
const C5 = zeros(5, 5)
const C6 = zeros(5, 5)
const C7 = zeros(5, 5)
const C8 = zeros(5, 5)
const C9 = zeros(5, 5)
const C10 = zeros(5, 5)
for i = 1:7
C4[i, i] = M_error[i] ^ 2
if i<=5
......@@ -58,9 +66,9 @@ for i = 1:7
C6[i, i] = ZS_over_ZP_err[i] ^ 2
C7[i, i] = ZV_err[i] ^ 2
C8[i, i] = RM_err[i] ^2
if i <= 4
C3[i, i] = t0_error[i] ^ 2
end
C9[i,i] = BA_MINUS_BP_err[i] ^2
C10[i,i] = ZERR[i] ^2
C3[i, i] = t0_error[i] ^ 2
end
end
......@@ -72,14 +80,18 @@ const a_ = t0_ph ./ sqrt.(8 .* t0_)
const M = cobs(M_values, C4, "charmed meson masses")
const Za = cobs(ZA_data, C5, "ZA")
const Zv = cobs(ZV_data, C7, "ZV")
const ZS_over_ZP = cobs(ZS_over_ZP_data, C6, "ZS over ZP")
const ZS_over_ZP = cobs(ZS_over_ZP_data, C6, "ZS/ZP")
const RM = cobs(RM_data, C8, "rm")
const BA_BP = cobs(BA_MINUS_BP_data, C9, "ba-bp")
const ZZ = cobs(ZDATA, C10, "Z")
zm(beta::Float64) = ZM[b_values .== beta][1]
zm_tm(beta::Float64) = ZM_tm[b_values .== beta][1]
t0(beta::Float64) = t0_[b_values2 .== beta][1]
a(beta::Float64) = a_[b_values2 .== beta][1]
t0(beta::Float64) = t0_[b_values .== beta][1]
a(beta::Float64) = a_[b_values .== beta][1]
za(beta::Float64) = Za[b_values .== beta][1]
zv(beta::Float64) = Zv[b_values .== beta][1]
zs_over_zp(beta::Float64) = ZS_over_ZP[b_values .== beta][1]
rm(beta::Float64) = RM[b_values .== beta][1]
\ No newline at end of file
rm(beta::Float64) = RM[b_values .== beta][1]
ba_bp(beta::Float64) = BA_BP[b_values .== beta][1]
Z(beta::Float64) = ZZ[b_values .== beta][1]
\ No newline at end of file
src/juobs.jl
View file @ f3ec800f
module juobs
using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim
using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim, ForwardDiff
using Base: Float64
import Statistics: mean
......
src/juobs_obs.jl
View file @ f3ec800f
......@@ -35,7 +35,7 @@ function meff(corr::Vector{uwreal}, plat::Vector{Int64}; pl::Bool=true, data::Bo
xlabel(L"$x_0$")
ylim(v-30*e, v+30*e)
# ylim(0.05, 0.075)
#xlim(64,length(y))
xlim(64,length(y))
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
......@@ -108,7 +108,7 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
# der_a0p = der_a0p + ca * der2_pp
#end
aux = der_a0p ./ ( corr_pp[3:end-2])
aux = der_a0p ./ ( 2. .* corr_pp[3:end-2])
mass = plat_av(aux, plat, wpm)
if pl
isnothing(wpm) ? uwerr(mass) : uwerr(mass, wpm)
......@@ -123,7 +123,7 @@ function mpcac(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}; ca:
errorbar(x, y, dy, fmt="x", color="black")
ylabel(L"$m_\mathrm{PCAC}$")
xlabel(L"$x_0$")
ylim(v-100*e, v+100*e)
ylim(v-20*e, v+20*e)
if !isnothing(kappa)
title(string(L"$\kappa_1 = $", kappa[1], L" $\kappa_2 = $", kappa[2]))
......@@ -330,13 +330,13 @@ function dec_const(vv::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64
corr_vv = vv[2:end-1]
T = length(corr_vv)
aux = exp.((collect(1:T) .- y0 ) .* fill(m, T))
aux = exp.((collect(1:T) .- y0 ) .* fill(m, T))
R = ((aux .* corr_vv).^2).^0.25
R_av = plat_av(R, plat, wpm)
f = sqrt(2 / m) * R_av
R .*= sqrt.(2 ./ [m for i in 1:length(R)])
if pl
R .*= sqrt.(2 ./ [m for i in 1:length(R)])
uwerr.(R)
R_av *= sqrt(2/m)
isnothing(wpm) ? uwerr(R_av) : uwerr(R_av, wpm)
......@@ -389,13 +389,13 @@ function dec_const_w(vv::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}
aux = exp.((collect(1:T-2) .- y0 ) .* fill(m, T-2))
else
#corr_vv = corr_vv[2:end-1]
aux = exp.((collect(1:T) .- y0 ) .* fill(m, T))
aux = exp.((collect(1:T) .- y0 ) .* fill(m/2, T))
end
R = (aux .* corr_vv)
R = (aux .* corr_vv ./ sqrt.(corr_pp))
R_av = plat_av(R, plat, wpm)
f = sqrt(2 / m) * sqrt(R_av)
f = sqrt(2 / m) * R_av
if pl
R .*= sqrt.(2 ./ [m for i in 1:length(R)])
uwerr.(R)
......
src/juobs_tools.jl
View file @ f3ec800f
......@@ -577,20 +577,23 @@ 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)
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, path_plt::Union{String,Nothing}=nothing,
plt_title::Union{Nothing,String}=nothing, label::Union{Nothing, LaTeXString}=nothing)
weight_model = Array{Float64,1}()
AIC = Array{Float64,1}()
chi2chi2exp = Array{Float64,1}()
p1 = Array{uwreal,1}()
mods = Array{String,1}()
Pval = Array{Float64,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) : [uwerr(y[i],wpm) for i in 1:length(y) ]
isnothing(wpm) ? uwerr.(y) : [uwerr(y[i],wpm) for i in 1:length(y)]
for INDEX in tmin_array ## vary tmin
for j in tmax_array ## vary tmax
......@@ -606,53 +609,112 @@ function bayesian_av(fun::Function, y::Array{uwreal}, tmin_array::Array{Int64},
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
chi2 = sum(fit.resid.^2) # * dof(fit) / chi_exp
Q = pvalue(chisq, sum(fit.resid.^2), value.(up), yy, wpm=wpm, W=W, nmc=10000)
# if chi2 / dof(fit) > 5. || abs(value(up[1])) > 2. || err(up[1])/value(up[1])*100 >= 10.
# error()
# end
push!(Pval, Q)
push!(AIC, chi2 + 2*k + 2*Ncut)
push!(chi2chi2exp, chi2 / dof(fit))
push!(p1, up[1])
push!(mods,string("[", INDEX+1, ",", j, "]"))
println("chi2_vs_exp " , chi2 / dof(fit), " chi2/dof ", sum(fit.resid.^2)/dof(fit))
# println("chi2_vs_exp " , chi2 / dof(fit), " chi2/dof ", sum(fit.resid.^2)/dof(fit))
catch e
@warn string(":/ Negative window for error propagation at tmin = ", INDEX, ", tmax = ", j, "; skipping that point")
end
end
end
# compute all weights
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)
# sort weights and discard 5% tail
idxW = sortperm(weight_model, rev=true)
cumulative_w = cumsum(weight_model[idxW])
# println("acceptance in bayeasian_av set to 0.99, should be 0.95")
idxcumw = findfirst(x->x>=0.95, cumulative_w)
# println(length(idxW) - idxcumw)
idxW = sort(idxW[1:idxcumw])
# take only the accepted
AIC = AIC[idxW]
Pval = Pval[idxW]
chi2chi2exp = chi2chi2exp[idxW]
mods = mods[idxW]
p1 = p1[idxW]
weight_model = weight_model[idxW]
weight_model ./= sum(weight_model)
p1_mean = sum(p1 .* weight_model)
isnothing(wpm) ? uwerr(p1_mean) : uwerr(p1_mean, wpm)
systematic_err = sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2)
if pl
if pl
fig = figure(figsize=(10,7.5))
subplots_adjust(hspace=0.1)
subplot(411)
if !isnothing(plt_title)
title(plt_title)
end
ax1 = gca()
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")
fill_between(x, v-e, v+e, color="royalblue", alpha=0.2)
errorbar(x, y, dy, fmt="^", mfc="none", color="navy", capsize=2)
setp(ax1.get_xticklabels(),visible=false) # Disable x tick labels
isnothing(label) ? ylabel(L"$p_1$") : ylabel(label)
display(gcf())
subplot(412)
ax2=gca()
bar(x, weight_model, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5)
setp(ax2.get_xticklabels(),visible=false) # Disable x tick labels
# errorbar(mods, weight_model, 0*dy, color="green")
ylabel(L"$W$")
subplot(413)
ax3=gca()
setp(ax3.get_xticklabels(),visible=false) # Disable x tick labels
bar(x, Pval, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5 )
ylabel(L"$p\ value$")
subplot(414)
bar(x, chi2chi2exp, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5 )
# ylabel(L"$\chi^2/\chi^2_{\mathrm{exp}}$")
ylabel(L"$\chi^2/dof$")
xticks(x, mods, rotation=45)
xlabel(L"$Models$")
display(fig)
if !isnothing(path_plt)
#tt = plt_title *".pdf"
savefig(path_plt)
end
close()
end
if !data
return (p1_mean, systematic_err)
else
return (p1_mean, systematic_err, p1, weight_model, AIC .+ offset)
FitStorage = Dict(
"W" => weight_model,
"AIC" => AIC .+ offset,
"pval" => Pval,
"chivsexp" => chi2chi2exp,
"mods" => mods
)
return (p1_mean, systematic_err, p1, FitStorage)
end
end
......@@ -702,7 +764,7 @@ function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_arra
push!(AIC, chi2 + 2*k2 + 2*Ncut)
push!(chi2chi2exp, chi2 / dof(fit))
push!(p1, up[1])
println(up[end])
push!(mods,string("[", INDEX+1, ",", j, "]"))
catch e
......@@ -750,20 +812,23 @@ function bayesian_av(fun1::Function, fun2::Function, y::Array{uwreal}, tmin_arra
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)
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,
path_plt::Union{String,Nothing}=nothing, plt_title::Union{Nothing,String}=nothing)
weight_model = Array{Float64,1}()
AIC = Array{Float64,1}()
chi2chi2exp = Array{Float64,1}()
p1 = Array{uwreal,1}()
mods = Array{String,1}()
AIC = Array{Float64,1}()
chi2chi2exp = Array{Float64,1}()
p1 = Array{uwreal,1}()
mods = Array{String,1}()
Pval = Array{Float64,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
isnothing(wpm) ? uwerr.(y) : [uwerr(y[i],wpm) for i in 1:length(y)]
for INDEX in tmin_array ## vary tmin
for j in tmax_array ## vary tmax
......@@ -781,10 +846,17 @@ function bayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{I
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, "]"))
Q = pvalue(chisq, sum(fit.resid.^2), value.(up), yy, wpm=wpm, W=W, nmc=10000)
if chi2 / dof(fit) < 5.
# if err(up[1])/value(up[1])*100 >= 10.
# error()
# end
push!(Pval, Q)
push!(AIC, chi2 + 2*k[indice] + 2*Ncut)
push!(chi2chi2exp, chi2 / dof(fit))
push!(p1, up[1])
push!(mods,string("[", INDEX+1, ",", j, "]"))
end
end
catch e
......@@ -792,42 +864,91 @@ function bayesian_av(fun::Array{Function}, y::Array{uwreal}, tmin_array::Array{I
end
end
end
# compute all weights
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)
# sort weights and discard 5% tail
idxW = sortperm(weight_model, rev=true)
cumulative_w = cumsum(weight_model[idxW])
idxcumw = findfirst(x->x>=0.95, cumulative_w)
println(length(idxW) - idxcumw)
idxW = sort(idxW[1:idxcumw])
# take only the accepted
AIC = AIC[idxW]
Pval = Pval[idxW]
chi2chi2exp = chi2chi2exp[idxW]
mods = mods[idxW]
p1 = p1[idxW]
weight_model = weight_model[idxW]
weight_model /= sum(weight_model)
p1_mean = sum(p1 .* weight_model)
isnothing(wpm) ? uwerr(p1_mean) : uwerr(p1_mean,wpm)
systematic_err = sqrt(sum(p1 .^ 2 .* weight_model) - (sum(p1 .* weight_model)) ^ 2)
if pl
if pl
subplots_adjust(hspace=0.1)
subplot(411)
if !isnothing(plt_title)
title(plt_title)
end
ax1 = gca()
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")
fill_between(x, v-e, v+e, color="royalblue", alpha=0.2)
errorbar(x, y, dy, fmt="^", mfc="none", color="navy", capsize=2)
setp(ax1.get_xticklabels(),visible=false) # Disable x tick labels
ylabel(L"$p_1$")
xlabel(L"model")
display(gcf())
subplot(412)
ax2=gca()
bar(x, weight_model, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5)
setp(ax2.get_xticklabels(),visible=false) # Disable x tick labels
# errorbar(mods, weight_model, 0*dy, color="green")
ylabel(L"$W$")
subplot(413)
ax3=gca()
setp(ax3.get_xticklabels(),visible=false) # Disable x tick labels
bar(x, Pval, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5 )
ylabel(L"$p\ value$")
subplot(414)
bar(x, chi2chi2exp, alpha=0.4, color="forestgreen", edgecolor="darkgreen", linewidth=1.5 )
ylabel(L"$\chi^2/\chi^2_{\mathrm{exp}}$")
xticks(x, mods, rotation=45)
xlabel(L"$Models$")
figure()
errorbar(mods, weight_model, 0*dy, color="green")
ylabel(L"$weight$")
xlabel(L"model")
display(gcf())
if !isnothing(path_plt)
tt = plt_title *".pdf"
savefig(joinpath(path_plt, tt))
end
close()
end
if !data
return (p1_mean, systematic_err)
else
return (p1_mean, systematic_err, p1, weight_model)
FitStorage = Dict(
"W" => weight_model,
"AIC" => AIC .+ offset,
"pval" => Pval,
"chivsexp" => chi2chi2exp,
"mods" => mods
)
return (p1_mean, systematic_err, p1, FitStorage)
end
end
......@@ -1296,4 +1417,123 @@ function inv_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Fl
cov = get_covar_multi_id(obs, wpm=wpm)
cov_inv = inv(cov)
return (cov_inv' + cov_inv) / 2.
end
@doc raw"""
pvalue(chisq::Function,
chi2::Float64,
xp::Vector{Float64},
data::Vector{uwreal},
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing} = Dict{Int64,Vector{Float64}}();
W::Union{Vector{Float64},Array{Float64,2}} = Vector{Float64}(),
nmc::Int64 = 5000)
Computes the p-value of a previously done fit, using as input the `\chi^2` observed from the fit, the fit parameters and the fitted data.
The p-value for a given `\chi^2` is the probability of, given the data you have, finding such a `\chi^2` or worse from a fit, and still
have the data well described by the fit function. `nmc` is the number of MC samples used to estimate the p-value integral, default is 5000.
By now it only works with a vector for weights (containing the diagonal of W)
```@example
function fit_defs(f::Function,x,W)
chisq(p,d)=sum((d-f(x,p)).^2 .*W)
return chisq
end
@.fun(x,p) = p[1] + p[2] * x
chisq = fit_defs(fun, x, 1.0 ./ err.(y) .^ 2)
fit = curve_fit(fun, x, value.(y), 1.0 ./ err.(y) .^ 2, [0.5, 0.5])
(up, chiexp) = fit_error(chisq, coef(fit), y)
wpm = Dict{Int64,Vector{Float64}}()
wpm[1] = [-1.0,-1.0,4-0,-1.0]
Q = pvalue(chisq, chi2, value.(up), y, wpm; W = 1.0 ./ err.(y) .^ 2, nmc=10000)
#Q = pvalue(chisq, chi2, value.(up), y; W = 1.0 ./ err.(y) .^ 2, nmc=10000)
#Q = pvalue(chisq, chi2, value.(up), y)
```
"""
function pvalue(chisq::Function,
chi2::Float64,
xp::Vector{Float64},
data::Vector{uwreal};
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing} = Dict{Int64,Vector{Float64}}(),
W::Union{Vector{Float64},Array{Float64,2}} = Vector{Float64}(),
nmc::Int64 = 5000)
n = length(xp) # Number of fit parameters
m = length(data) # Number of data
xav = zeros(Float64, n+m)
for i in 1:n
xav[i] = xp[i]
end
for i in n+1:n+m
xav[i] = data[i-n].mean
end
ccsq(x::Vector) = chisq(view(x, 1:n), view(x, n+1:n+m))
if (n+m < 4)
cfg = ForwardDiff.HessianConfig(ccsq, xav, ADerrors.Chunk{1}());
else
cfg = ForwardDiff.HessianConfig(ccsq, xav, ADerrors.Chunk{4}());
end
hess = Array{Float64}(undef, n+m, n+m)
ForwardDiff.hessian!(hess, ccsq, xav, cfg)
cse = 0.0
Q = dQ = 0.0
if (m-n > 0)
if (length(W) == 0)
Ww = zeros(Float64, m)
for i in 1:m
if (data[i].err == 0.0)
#isnothing(wpm) ? wuerr(data[i]) : uwerr(data[i], wpm)
uwerr(data[i], wpm)
if (data[i].err == 0.0)
error("Zero error in fit data")
end
end
Ww[i] = 1.0 / data[i].err^2
end
else
Ww = W
end
#cse = chiexp(hess, data, Ww, wpm)
m = length(data)
n = size(hess, 1) - m
hm = view(hess, 1:n, n+1:n+m)
sm = Array{Float64, 2}(undef, n, m)
for i in 1:n, j in 1:m
sm[i,j] = hm[i,j] / sqrt.(Ww[j])
end
maux = sm * sm'
hi = LinearAlgebra.pinv(maux)
Px = -hm' * hi * hm
for i in 1:m
Px[i,i] = Ww[i] + Px[i,i]
end
C = cov(data)
nu = sqrt(C) * Px * sqrt(C)
N = length(nu[1,:])
z = randn(N, nmc)
eig = abs.(eigvals(nu))
eps = 1e-14 * maximum(eig)
eig = eig .* (eig .> eps)
aux = eig' * (z .^ 2)
Q = 1.0 - juobs.mean(aux .< chi2)
x = chi2 .- eig[2:end]' * (z[2:end,:].^2)
x = x / eig[1]
#dQ = juobs.mean((x .> 0) .* exp.(-x * 0.5) * 0.5 ./ sqrt.(abs.(x)))
#dQ = err(cse)/value(cse) * dQ
end
return Q #uwreal([Q,dQ],"")
end
\ No newline at end of file
src/juobs_types.jl
View file @ f3ec800f
......@@ -182,7 +182,7 @@ mutable struct Corr
theta2 = h[1].theta2
return new(a, kappa, mu, gamma, y0, theta1, theta2)
end
Corr(o::Vector{uwreal}, k::Vector{Float64}, m::Vector{Float64}, g::Vector{String}, y::Int64) = new(o, k, m, g, y)
Corr(o::Vector{uwreal}, k::Vector{Float64}, m::Vector{Float64}, g::Vector{String}, y::Int64, th1::Vector{Float64}, th2::Vector{Float64} ) = new(o, k, m, g, y, th1, th2)
end
mutable struct YData
......
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