pvalue now has a keyword option C::AbstractMatrix{Float64} defaulted to ADerrors.cov(data,wmp)

src/juobs_tools.jl

@@ -1235,7 +1235,7 @@ function get_chi(model::Function,x::Vector,par,data,W::Vector{Float64})
 end
 
 @doc raw"""
-    fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, npar::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),corr::Bool = true, W::VecOrMat{Float64}=Vector{Float64}(),guess::Vector{Float64} = fill(0.5,npar));
+    fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, npar::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}}=Dict{Int64,Vector{Float64}}(),corr::Bool = true, W::VecOrMat{Float64}=Vector{Float64}(),guess::Vector{Float64} = fill(0.5,npar),logfile);
 
 It executes a fit of the `ydata` with the fit function `model`.
 
@@ -1246,11 +1246,13 @@ It executes a fit of the `ydata` with the fit function `model`.
   - `W::VecOrMat{Float64}`: Weight matrix/vector. When given, the flag `corr` is overridden. If `W` is a `Vector` the fit is perfomed uncorrelated,otherwise is correlated. 
   - `guess::Vector{Float64}`: Initial guess used in the fit routine.
   - `logfile`: handle to a logfile. If `nothing`, no log will be written. It can be `IO` file or a custom log structure that has an overloading for `print()` and `println()`.
+
+
 # Returns
 It returns a NamedTuple with names:
   - `:par`:    fit parameter as `Vector{uwreal}`
   - `:chi2`:   chisquare,
-  - `:chiexp`: chisquare expected
+  - `:chiexp`: chisquare expected 
   - `:pval`:   pvalue 
 """
 function fit_routine(model::Function,
@@ -1261,7 +1263,7 @@ function fit_routine(model::Function,
   corr::Bool = true, 
   W::VecOrMat{Float64}=Vector{Float64}(),
   guess::Vector{Float64} = fill(0.5,npar),
-  logfile = nothing);
+  logfile = nothing,);
 
   nobs = length(ydata)
   if length(xdata)!=nobs
@@ -1271,7 +1273,7 @@ function fit_routine(model::Function,
   [uwerr(y, wpm) for y in ydata]
   if length(W) ==0
     if corr 
-      W = ADerrors.cov(ydata,wpm)
+      C = ADerrors.cov(ydata,wpm)
       W = LinearAlgebra.pinv(W);
       W = (W'+W)*0.5
     else
@@ -1283,21 +1285,26 @@ function fit_routine(model::Function,
 
   chisq(p,d) = get_chi(model,xdata,p,d,W)
   chi2 = sum(fit.resid.^2)
-  par,chiexp = fit_error(chisq,LsqFit.coef(fit),ydata,wpm,W)
+  par = fit_error(chisq,LsqFit.coef(fit),ydata,wpm,W,false)
+
+  chiexp = ADerrors.chiexp(chisq,LsqFit.coef(fit),ydata,wpm,W=W)
+
   pval = juobs.pvalue(chisq,chi2,LsqFit.coef(fit),ydata,W)
 
-  if !isnothing(logfile) 
-    println(logfile, "Fit routine in [$(xdata[1]),...,$(xdata[2])]")
-    println(logifle, "parameters :", par...)
-    println(logfile, "chi2:       ", chi2)
-    println(logfile, "chiexp:     ", chiexp)
-    println(logfile, "pvalue:     ", pval)
+  if !isnothing(logfile)
+
+    println(logfile, "juobs.fit_routine output:") 
+    println(logfile, "\t\tFit routine in [$(xdata[1]),...,$(xdata[2])]")
+    println(logfile, "\t\tparameters :", par...)
+    println(logfile, "\t\tchi2:       ", chi2)
+    println(logfile, "\t\tchiexp:     ", chiexp)
+    println(logfile, "\t\tpvalue:     ", pval)
   end
 
   return (par=par,chi2=chi2,chiexp=chiexp,pval=pval)
 end
 
-function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
+#= function fit_routine(model::Vector{Function}, xdata::Vector{Array{Float64, N}} where N, ydata::Vector{Array{uwreal, N}} where N, param::Int64; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing,
     correlated_fit::Bool=false)
 
     if !(length(model) == length(xdata) == length(ydata))
@@ -1800,7 +1807,7 @@ function inv_covar_multi_id(obs::Vector{uwreal}; wpm::Union{Dict{Int64,Vector{Fl
     cov = get_covar_multi_id(obs, wpm=wpm)
     cov_inv = inv(cov)
     return (cov_inv' + cov_inv) / 2.
-end
+end =#
 
 @doc raw"""
     pvalue(chisq::Function,
@@ -1840,7 +1847,8 @@ function pvalue(chisq::Function,
     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)
+    nmc::Int64 = 5000,
+    C::AbstractMatrix{Float64} = ADerrors.cov(data,wpm))
 
     n = length(xp)   # Number of fit parameters
     m = length(data) # Number of data
@@ -1892,8 +1900,6 @@ function pvalue(chisq::Function,
         Px[i,i] = W[i] + Px[i,i]
     end
 
-    C = cov(data) 
-        
     nu = sqrt(C) * Px * sqrt(C)
         
     N = length(nu[1,:])
@@ -1919,7 +1925,8 @@ function pvalue(chisq::Function,
     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)
+    nmc::Int64 = 5000,
+    C::AbstractMatrix{Float64} = ADerrors.cov(data,wpm))
 
     n = length(xp)   # Number of fit parameters
     m = length(data) # Number of data
@@ -1955,7 +1962,7 @@ function pvalue(chisq::Function,
         hi   = LinearAlgebra.pinv(maux)
         Px   = W - hm' * hi * hm
         
-        C = cov(data) 
+        
         
         nu = sqrt(C) * Px * sqrt(C)