bug in pval, updated to support non-diagonal matrices as weights

src/juobs_tools.jl

@@ -1588,6 +1588,153 @@ Q = pvalue(chisq, chi2, value.(up), y, wpm; W = 1.0 ./ err.(y) .^ 2, nmc=10000)
 #Q = pvalue(chisq, chi2, value.(up), y)
 ```
 """
+
+function juobs.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::Vector{Float64} = 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)
+
+    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
+            W = Ww
+        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])
+        end
+        maux = sm * sm'
+        hi   = LinearAlgebra.pinv(maux)
+        Px   = -hm' * hi * hm
+
+        for i in 1:m
+            Px[i,i] = W[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
+end
+
+function juobs.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::Array{Float64,2},
+    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)
+
+    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
+
 function pvalue(chisq::Function,
     chi2::Float64,
     xp::Vector{Float64}, 
@@ -1642,8 +1789,14 @@ function pvalue(chisq::Function,
         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])
+        if typeof(Ww) == Array{Float64,2}
+            for i in 1:n, j in 1:m
+                sm[i,j] = hm[i,j] .* sqrt(inv(Ww))[i,j]
+            end
+        elseif typeof(Ww) == Vector{Float64}
+            for i in 1:n, j in 1:m
+                sm[i,j] = hm[i,j] / sqrt.(Ww[j])
+            end
         end
         maux = sm * sm'
         hi   = LinearAlgebra.pinv(maux)