Commit 1780005b authored by Antonino D'Anna's avatar Antonino D'Anna

updated documentation on pvalue. Slight changes to make it more intuitive to call pvalue

parent 1c430848
......@@ -183,7 +183,7 @@ function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bo
if a0p.kappa == pp.kappa || a0p.mu == pp.mu
if a0p.mu == [0.0, 0.0]
return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm,savepl=davepl, label=label)
return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, kappa=a0p.kappa, mu=nothing, wpm=wpm,savepl=savepl, label=label)
else
return mpcac(a0p.obs, pp.obs, plat, ca=ca, pl=pl, data=data, mu=a0p.mu, kappa=nothing, wpm=wpm,savepl=savepl,label=label)
end
......@@ -191,7 +191,7 @@ function mpcac(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0, pl::Bo
error("mu or kappa values does not match")
end
end
#=
@doc raw"""
mpcac_ren(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, za_zp::uwreal; ca::Float64=0.0,ba_bp_z::uwreal=0.0,pl::Bool=true, data::Bool=false,
kappa::Union{Vector{Float64}, Nothing}=nothing, mu::Union{Vector{Float64}, Nothing}=nothing,
......@@ -313,7 +313,7 @@ function mpcac_ren(a0p::Corr, pp::Corr, plat::Vector{Int64}; ca::Float64=0.0,ba_
error("mu or kappa values does not match")
end
end
=#
## Decay constants
@doc raw"""
dec_const(a0p::Vector{uwreal}, pp::Vector{uwreal}, plat::Vector{Int64}, m::uwreal, y0::Int64; ca::Float64=0.0, pl::Bool=true, data::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
......
......@@ -1490,18 +1490,25 @@ function inv_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Fl
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)
function pvalue(chisq::Function,
chi2::Float64,
xp::Vector{Float64},
data::Vector{uwreal},
W::Union{Vector{Float64},Nothing} = nothing,
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing} = Dict{Int64,Vector{Float64}}();
nmc::Int64 = 5000)
function pvalue(chisq::Function,
chi2::Float64,
xp::Vector{Float64},
data::Vector{uwreal},
W::Array{Float64,2};
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing} = Dict{Int64,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)
......@@ -1517,16 +1524,16 @@ fit = curve_fit(fun, x, value.(y), 1.0 ./ err.(y) .^ 2, [0.5, 0.5])
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, va(Nx − NA )/2(Nx − NA )/2lue.(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}}(),
data::Vector{uwreal},
W::Union{Vector{Float64},Nothing} = nothing;
wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing} = nothing,
nmc::Int64 = 5000)
n = length(xp) # Number of fit parameters
......@@ -1550,18 +1557,19 @@ function pvalue(chisq::Function,
ForwardDiff.hessian!(hess, ccsq, xav, cfg)
if (m-n > 0)
W = zeros(Float64, m)
if !isnothing(W)
W = 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)
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
global W[i] = 1.0 / data[i].err^2
W[i] = 1.0 / data[i].err^2
end
end
m = length(data)
n = size(hess, 1) - m
......
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