pvalue function for uncorrelated data now accept a W vector. default value is...

pvalue function for uncorrelated data now accept a W vector. default value is Vector{Float64}(), and if W is not provided, it is computed as the inverse of the error squared of the data.

src/Changelog.md

@@ -9,4 +9,5 @@ This file should not be merged into master, but only kept as a reference within
     - Updated documentation in juobs_tools.jl
     - Updated documentation in juobs_type.jl 
   - Dic 8 2024: added functions plat_av to the existing one. This change should not be code-breaking, since it only add different version of the exising one. Now it is possible to pass a Vector or a Matrix W containing the weight to use in the plateau average. If W is a Matrix, a correlated fit is perform. If W is not given, then the old function is called. 
-  - Dic 10 2024: mpcac function now follow the convention more clearly. documentation updated juobs_tools
+  - Dic 10 2024: mpcac function now follow the convention more clearly. documentation updated juobs_tools. 
\ No newline at end of file
+                 pvalue function for uncorrelated data now accept a W vector. default value is Vector{Float64}(), and if W is not provided, it is computed as the inverse of the error squared of the data. 
\ No newline at end of file

src/juobs_tools.jl

@@ -1838,12 +1838,16 @@ Q = pvalue(chisq, chi2, value.(up), y, wpm; W = 1.0 ./ err.(y) .^ 2, nmc=10000)
 function pvalue(chisq::Function,
     chi2::Float64,
     xp::Vector{Float64}, 
-    data::Vector{uwreal};
+    data::Vector{uwreal},
+    W::Vector{Float64}=Vector{Float64}();
     wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing} =  Dict{Int64,Vector{Float64}}(),
     nmc::Int64 = 5000)
 
     n = length(xp)   # Number of fit parameters
     m = length(data) # Number of data
+    if (m-n<0)
+      error("more parameters than data")
+    end
 
     xav = zeros(Float64, n+m)
     for i in 1:n
@@ -1861,55 +1865,52 @@ function pvalue(chisq::Function,
         
     hess = Array{Float64}(undef, n+m, n+m)
     ForwardDiff.hessian!(hess, ccsq, xav, cfg)
-
-    if (m-n > 0)
-        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)
-                    if (data[i].err == 0.0)
-                        error("Zero error in fit data")
-                    end
-                end
-                global W[i] = 1.0 / data[i].err^2
-            end
-
-        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.(W[j])
+    
+    if length(W) ==0 
+      W = zeros(Float64, m)
+      for i in 1:m
+        if (data[i].err == 0.0)
+          uwerr(data[i], wpm)
+          if (data[i].err == 0.0)
+            error("Zero error in fit data")
+          end
         end
-        maux = sm * sm'
-        hi   = LinearAlgebra.pinv(maux)
-        Px   = -hm' * hi * hm
+        W[i] = 1.0 / data[i].err^2
+      end
+    end
 
-        for i in 1:m
-            Px[i,i] = W[i] + Px[i,i]
-        end
+    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.(W[j])
+    end
+    maux = sm * sm'
+    hi   = LinearAlgebra.pinv(maux)
+    Px   = -hm' * hi * hm
 
-        C = cov(data) 
+    for i in 1:m
+        Px[i,i] = W[i] + Px[i,i]
+    end
+
+    C = cov(data) 
         
-        nu = sqrt(C) * Px * sqrt(C)
+    nu = sqrt(C) * Px * sqrt(C)
         
-        N = length(nu[1,:])
-        z = randn(N, nmc)
+    N = length(nu[1,:])
+    z = randn(N, nmc)
 
-        eig = abs.(eigvals(nu))
-        eps = 1e-14 * maximum(eig)
-        eig = eig .* (eig .> eps)
+    eig = abs.(eigvals(nu))
+    eps = 1e-14 * maximum(eig)
+    eig = eig .* (eig .> eps)
 
-        aux = eig' * (z .^ 2)
-        Q = 1.0 - juobs.mean(aux .< chi2)
+    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
+    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
     return Q
 end