Including Mrat covariances in juobs_const.

Study and preliminary tests  on correlated fits

constants/juobs_const.jl

@@ -24,6 +24,14 @@ const t0_error = [11, 16, 18, 29] .* 1e-3
 const t0_ph_value = [0.413]
 const t0_ph_error = ones(1,1) .* 5e-3
 
+## taking into account correlations in zm_tm
+
+C = [[0.375369e-6, 0.429197e-6, -0.186896e-5] [0.429197e-6, 0.268393e-4, 0.686776e-4] [-0.186896e-5, 0.686776e-4, 0.212386e-3]]
+z = cobs([0.348629, 0.020921, 0.070613], C, "z")
+ZP(b) = z[1] + z[2] * (b - 3.79) + z[3] * (b - 3.79)^2
+Mrat = uwreal([.9148, 0.0088], "M/mhad")
+zm_tm_covar(b) = Mrat / ZP(b)
+##
 #Cov matrices
 const C1 = zeros(5, 5)
 const C2 = zeros(5, 5)

src/juobs.jl

@@ -10,7 +10,7 @@ include("juobs_tools.jl")
 include("juobs_obs.jl")
 
 export read_mesons, read_ms1, read_ms, read_md, truncate_data!
-export get_matrix, energies, uwdot, uweigvals, uweigvecs, uweigen, invert, getall_eigvals, getall_eigvecs, hess_reduce, uwcholesky, transpose, tridiag_reduction
+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_sym, md_sea, md_val, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine, global_fit_routine
 export meff, mpcac, dec_const, dec_const_pcvc, comp_t0
 

src/juobs_linalg.jl

@@ -741,6 +741,33 @@ function idty(n::Int64)
     return res
 end
 
+function make_positive_def(A::Matrix)
+    B = 0.5 * (A + A')
+    U, S, V = svd(B)
+    H = V * Diagonal(S) * V'
+    A_hat = 0.5 * (B + H)
+    A_hat = 0.5 * ( A_hat + A_hat')
+    k = 0
+    while !isposdef(A_hat)
+        mineig = minimum(eigvals(A_hat))
+        eps =  2.220446049250313e-16
+        A_hat = A_hat + (-mineig*k^2 .+ eps*Matrix(I, size(A)))
+        k+=1
+    end
+    return A_hat
+end
+
+function invert_covar_matrix(A::Matrix)
+    F = svd(A)
+    cov_inv = F.V * Diagonal(1 ./F.S) * F.U'
+    return cov_inv  
+end
+
+function more_penrose_inv(A::Matrix)
+    F = svd(A)
+
+end
+
 
 
 """

src/juobs_obs.jl

@@ -569,7 +569,7 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false,
     
     model(x, p) = get_model(x, p, npol)
 
-    par = fit_routine(model, x, t2E, npol)
+    par, chi  = fit_routine(model, x, t2E, npol)
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl
@@ -661,7 +661,7 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
     
     model(x, p) = get_model(x, p, npol)
 
-    par = fit_routine(model, x, t2E, npol)
+    par, chi = fit_routine(model, x, t2E, npol)
     fmin(x, p) = model(x, p) .- 0.3
     t0 = root_error(fmin, t[nt0], par)
     if pl

src/juobs_tools.jl

@@ -382,25 +382,44 @@ fit_routine(model, xdata, ydata, param=3)
 fit_routine(model, xdata, ydata, param=3, covar=true)
 ```
 """
-function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; info::Bool=false, wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing, correlated_fit::Bool=false)
     isnothing(wpm) ? uwerr.(ydata) : [uwerr(yaux, wpm) for yaux in ydata]
     
     yval = value.(ydata)
     yer = err.(ydata)
-    # Generate chi2 + solver
+    if !correlated_fit
         chisq = gen_chisq(model, xdata, yer)
         fit = curve_fit(model, xdata, yval, 1.0 ./ yer.^2, fill(0.5, param))
         (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata) : fit_error(chisq, coef(fit), ydata, wpm)
+    else
+        covar_inv = make_positive_def(invert_covar_matrix(cov(ydata))) 
+        #covar_inv = (covar_inv + covar_inv')/2
+        chisq(par,dat) = gen_chisq_correlated(model, xdata, covar_inv, par, dat)
+        fit = curve_fit(model, xdata, yval, covar_inv, fill(0.5, param))        
+        println(chisq(coef(fit), ydata))
+        (upar, chi_exp) = isnothing(wpm) ? fit_error(chisq, coef(fit), ydata, W=covar_inv) : fit_error(chisq, coef(fit), ydata, wpm, W=covar_inv)
+
+    end
+    chi2_fit_res = sum(fit.resid.^2 )
+    println("chi2 from fit residual = ", chi2_fit_res)
+    println("chi2 from chi2 function = ", chisq(coef(fit), ydata))
+
+
     #Info
+    if info
         for i = 1:length(upar)
             isnothing(wpm) ? uwerr(upar[i]) : uwerr(upar[i], wpm)
             print("\n Fit parameter: ", i, ": ")
             details(upar[i])
         end
-    chis2_corrected = (length(yval) - param) * chisq(coef(fit), ydata) / chi_exp
-    println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
+    end
+    #chis2_corrected = (length(yval) - param) * chisq(coef(fit), ydata) / chi_exp
+    chis2_corrected = (length(yval) - param) * chi2_fit_res / chi_exp
+
+    #println("Chisq / chiexp: ", chisq(coef(fit), ydata), " / ", chi_exp, " (dof: ", length(yval) - param,")")
+    println("Chisq / chiexp: ", chi2_fit_res, " / ", chi_exp, " (dof: ", length(yval) - param,")")
     println("Chisq corrected: ", chis2_corrected)
-    return upar, value(chis2_corrected)
+    return upar, chi2_fit_res / chi_exp
 end
 
 function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; 
@@ -457,6 +476,8 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         println("Chisq / chiexp: ", min_fun_cov(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
         chis2_corrected = (length(ydata) - param) * min_fun_cov(sol.minimizer) / chi2_exp
         println("Chisq corrected: ", chis2_corrected)
+        return upar, min_fun_cov(sol.minimizer) / chi2_exp
+
     else
         chisq_full(p, d) = get_chi2(model, d, ddat, p, Nalpha)
         min_fun(t) = chisq_full(t, dat)
@@ -466,6 +487,7 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
         println("Chisq / chiexp: ", min_fun(sol.minimizer), " / ", chi2_exp, " (dof: ", length(ydata) - param,")")
         chis2_corrected = (length(ydata) - param) * min_fun(sol.minimizer) / chi2_exp
         println("Chisq corrected: ", chis2_corrected)
+        return upar, min_fun(sol.minimizer) / chi2_exp
     end
 
     #### chisq_full, min_fun out of conditional ->
@@ -477,7 +499,6 @@ function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}
     #    print("\n Fit parameter: ", i, ": ")
     #    details(upar[i])
     #end
-    return upar, chis2_corrected
     
 end
 
@@ -505,7 +526,6 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6
         fx = similar(x_flat)
         for ii in 1:N_sets
             mod = models[ii]
-            #par = param[idx[ii]]             
             fx[chunk_idx[ii]] = mod(x_flat[chunk_idx[ii]], param)
         end
         return fx
@@ -514,26 +534,13 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6
         chi = 0.0
         for ii in 1:N_sets
             mod = models[ii]
-            #par = param[idx[ii]]
             # uncorrelated chi2
             chi += sum((data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param)).^2 ./ err.(y_flat[chunk_idx[ii]]).^2 )        
             # correlated chi2
-            #cov_mat = cov(data[chunk_idx[ii]])
-            #chi += sum(
-            #    juobs.transpose(data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param)) * cov_mat^(-1) * (data[chunk_idx[ii]] .- mod(x_flat[chunk_idx[ii]], param))
-            #)
         end
         return chi
     end
-    #tmp_par = [0.1, 0.6, 1, 1, 0.5, 0.75, 1, 1, 1, 1, 1, 1]
-    #tmp_par = [1, 0.1, 1, 0.6]
-    #tmp_par = [1, 0.5, 1, 0.75]
-    #tmp = flat_chisq(tmp_par, y_flat)
-    #println("chi2 with correct values = ", tmp)
     fit = curve_fit(flat_model, x_flat, value.(y_flat), 1.0 ./ err.(y_flat).^2, guess_param )
-    #W = cov(y_flat)^(-1)
-    #println(typeof(W))
-    #fit = curve_fit(flat_model, x_flat, value.(y_flat), W, guess_param )
     #uncorrelated chi2exp
     (upar, chi2_exp) = isnothing(wpm) ? fit_error(flat_chisq, coef(fit) , y_flat) : fit_error(flat_chisq, coef(fit) , y_flat, wpm)
     #correlated chi2exp
@@ -541,7 +548,6 @@ function global_fit_routine(models::Vector{Function}, xdata::Vector{Array{Float6
     chisq = value(flat_chisq(coef(fit), y_flat))  
     if log
         println("Converged param: ", fit.converged)
-        #println("wt param: ", fit.wt)
         println("Chisq / chiexp: ", chisq, " / ", chi2_exp, " (dof: ", dof(fit),")")
         chis2_corrected = (dof(fit)) * chisq / chi2_exp
         println("Chisq corrected: ", chis2_corrected)
@@ -553,6 +559,24 @@ function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64}) #constra
     chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
     return chisq
 end
+function gen_chisq_correlated(f::Function, x::Array{<:Real}, covar_inv::Matrix, par, dat )
+    chi2 = 0.0
+    #delta = dat - f(x, par)
+    #chi2 = (reshape(delta, 1,:) * covar_inv * delta )[1]
+    for i in 1:length(dat)
+        for j in 1:length(dat)
+            #println(f(x,par)[i])
+            #println(f(x,par)[j])
+            #println(f(x[j,:],par))
+            #println(f(x[i,:],par))
+            chi2 = chi2 +  (dat[i] - f(x,par)[i]) * covar_inv[i,j] * (dat[j] - f(x,par)[j])
+        end
+    end
+    #chisq(par, dat) = sum(
+    #    transpose((dat .- f(x,par))) * covar_inv * (dat .- f(x,par))
+    #) 
+    return chi2
+end
 
 function get_chi2(f::Function, data, ddata, par, Nalpha) #full
     chi2 = 0.0