Added pvalue from Alejandro while waiting for merge

src/juobs_tools.jl

@@ -1521,12 +1521,12 @@ Q = pvalue(chisq, chi2, value.(up), y, wpm; 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
@@ -1549,11 +1549,8 @@ function pvalue(chisq::Function,
     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)
+        W = zeros(Float64, m)
             for i in 1:m
                 if (data[i].err == 0.0)
                     #isnothing(wpm) ? wuerr(data[i]) : uwerr(data[i], wpm)
@@ -1562,29 +1559,25 @@ function pvalue(chisq::Function,
                         error("Zero error in fit data")
                     end
                 end
-                Ww[i] = 1.0 / data[i].err^2
-            end
-        else
-            Ww = W
+                global W[i] = 1.0 / data[i].err^2
             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])
+            sm[i,j] = hm[i,j] / sqrt.(W[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]
+            Px[i,i] = W[i] + Px[i,i]
         end
+
         C = cov(data) 
         
         nu = sqrt(C) * Px * sqrt(C)
@@ -1604,6 +1597,69 @@ function pvalue(chisq::Function,
         #dQ = juobs.mean((x .> 0) .* exp.(-x * 0.5) * 0.5 ./ sqrt.(abs.(x)))
         #dQ = err(cse)/value(cse) * dQ
     end
+    return Q
+end
+
+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)
+
+    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
         
-    return Q #uwreal([Q,dQ],"")
+    hess = Array{Float64}(undef, n+m, n+m)
+    ForwardDiff.hessian!(hess, ccsq, xav, cfg)
+
+    if (m-n > 0)
+        m = length(data)
+        n = size(hess, 1) - m
+        
+        Lm = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(W))
+        Li = LinearAlgebra.inv(Lm.L)
+        
+        hm = view(hess, 1:n, n+1:n+m)
+        sm = hm * Li'
+        
+        maux = sm * sm'
+        hi   = LinearAlgebra.pinv(maux)
+        Px   = W - hm' * hi * hm
+        
+        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
 end