pvalue upd

src/juobs.jl

@@ -11,7 +11,7 @@ include("juobs_obs.jl")
 
 export read_mesons, read_mesons_correction, read_ms1, read_ms, read_md, truncate_data!, concat_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, bayesian_av
+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, pvalue
 export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 
 end # module

src/juobs_tools.jl

@@ -1187,21 +1187,39 @@ end
     chi2::Float64,
     xp::Vector{Float64}, 
     data::Vector{uwreal},
-    wpm::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}();
+    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.
+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)
 
-Q = pvalue(chisq, chi2, value.(up), y, wpm; W=W, nmc=10000)
 """
 function pvalue(chisq::Function,
     chi2::Float64,
     xp::Vector{Float64}, 
     data::Vector{uwreal},
-    wpm::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}();
+    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)