Skip to content

J
juobs
Commit 5f3ab2cc authored by Antonino's avatar Antonino
Browse files
Options

Reverted Gamma method to accept only Vector

parent 30cb3eb7
Hide whitespace changes
Inline Side-by-side
Showing with 171 additions and 171 deletions
src/juobs_gammaMethod.jl
View file @ 5f3ab2cc
"""
GammaMethod
It includes function for the the Gamma Method. It is intended to complement ADerrors, so any function that is made obsolete
by an newer version of ADerrors should be deprecrated.
It includes function for the the Gamma Method. It is intended to complement ADerrors, so any function that is made obsolete
by an newer version of ADerrors should be deprecrated.
# Methods
# Methods
- calc_gamma(a1::ADerrors.uwreal [,a2::Aderrors.uwreal] ,id::Union{Int64,String}[, ws::ADerrors.wspace])
- bias(a1::ADerrors.uwreal [,a2::ADerrors.uwreal], id::Union{Int64,String} [, ws::Aderrors.wspace]; iw::Int64=0)
- cov_pe(vobs::Vectir{ADerrors.uwreal}[,ws:ADerrors.wspace])
......@@ -14,11 +14,11 @@ by an newer version of ADerrors should be deprecrated.
- chiexp(chisq::Function,par,d,C,W)
"""
module GammaMethod
import ADerrors, ForwardDiff, LinearAlgebra
"""
import ADerrors, ForwardDiff, LinearAlgebra
"""
calc_gamma(a1::ADerrors.uwreal [,a2::Aderrors.uwreal] ,id::Union{Int64,String}[, ws::ADerrors.wspace])
Computes the Gamma function associated to the ensemble `id` without the bias included by ADerrors for t>=0
If `id` is not associated to a Montecarlo chain, it return a Vector with the correlation between `a1` and `a2`
......@@ -31,7 +31,7 @@ function calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerro
if !(id in a1.ids) || !(id in a2.ids)
return [0.0]
end
if nd==1
function get_v(a,id=id,ws=ws)
v = 0
......@@ -51,15 +51,15 @@ function calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerro
nrep = length(ws.fluc[idx].ivrep)
ftemp1 = Dict{Int64,Vector{Complex{Float64}}}()
ftemp2 = Dict{Int64,Vector{Complex{Float64}}}()
for i in 1:nrep
ns = length(ws.fluc[idx].fourier[i])
ftemp1[i] = zeros(Complex{Float64}, ns)
ftemp2[i] = zeros(Complex{Float64}, ns)
end
gamma = zeros(Float64, nt)
for i in eachindex(a1.prop)
if (a1.prop[i] && (ws.map_nob[i] == id))
for k in 1:nrep
......@@ -69,62 +69,62 @@ function calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerro
end
for i in eachindex(a2.prop)
if (a2.prop[i] && (ws.map_nob[i] == id))
for k in 1:nrep
ftemp2[k] = ftemp2[k] + a2.der[i]*ws.fluc[i].fourier[k]
if (a2.prop[i] && (ws.map_nob[i] == id))
for k in 1:nrep
ftemp2[k] = ftemp2[k] + a2.der[i]*ws.fluc[i].fourier[k]
end
end
end
end
for k in 1:nrep
ftemp2[k] .= ftemp1[k].*conj(ftemp2[k])
ADerrors.FFTW.ifft!(ftemp2[k])
for ig in 1:min(nt,length(ftemp2[k]))
gamma[ig] = gamma[ig] + real(ftemp2[k][ig])
end
ftemp2[k] .= ftemp1[k].*conj(ftemp2[k])
ADerrors.FFTW.ifft!(ftemp2[k])
for ig in 1:min(nt,length(ftemp2[k]))
gamma[ig] = gamma[ig] + real(ftemp2[k][ig])
end
end
for ig in 1:nt
nrcnt = count(map(x -> div(x, 2*ws.fluc[idx].ibn), ws.fluc[idx].ivrep) .> ig-1)
gamma[ig] = gamma[ig] / (nd_eff - nrcnt*(ig-1))
end
return gamma
end
calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String,ws::ADerrors.wspace) = calc_gamma(a1,a2,ws.str2id[id],ws)
calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String) = calc_gamma(a1,a2,ADerrors.wsg.str2id[id],ADerrors.wsg)
calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64) = calc_gamma(a1,a2,id,ADerrors.wsg)
calc_gamma(a::ADerrors.uwreal,id::Int64,ws::ADerrors.wspace) = calc_gamma(a,a,id,ws)
calc_gamma(a::ADerrors.uwreal,id::String,ws::ADerrors.wspace) = calc_gamma(a,a,ws.str2id[id],ws)
calc_gamma(a::ADerrors.uwreal,id::String) = calc_gamma(a,a,ADerrors.wsg.str2id[id],ADerrors.wsg)
calc_gamma(a::ADerrors.uwreal,id::Int64) = calc_gamma(a,a,id,ADerrors.wsg)
function drho(rho::Vector{Float64},nd::Int64,iw::Int64=0)
return gamma
end
calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String,ws::ADerrors.wspace) = calc_gamma(a1,a2,ws.str2id[id],ws)
calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String) = calc_gamma(a1,a2,ADerrors.wsg.str2id[id],ADerrors.wsg)
calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64) = calc_gamma(a1,a2,id,ADerrors.wsg)
calc_gamma(a::ADerrors.uwreal,id::Int64,ws::ADerrors.wspace) = calc_gamma(a,a,id,ws)
calc_gamma(a::ADerrors.uwreal,id::String,ws::ADerrors.wspace) = calc_gamma(a,a,ws.str2id[id],ws)
calc_gamma(a::ADerrors.uwreal,id::String) = calc_gamma(a,a,ADerrors.wsg.str2id[id],ADerrors.wsg)
calc_gamma(a::ADerrors.uwreal,id::Int64) = calc_gamma(a,a,id,ADerrors.wsg)
function drho(rho::Vector{Float64},nd::Int64,iw::Int64=0)
iw = iw ==0 ? ADerrors.wopt_ulli(nd,ADerrors.DEFAULT_STAU,rho) : iw
nt = length(rho);
dr = zeros(nt);
for t in 1:nt
is = max(1, t-iw-2) + 1
ie = is + iw-1
for m in is:ie
aux=0.0
if m <= nt
aux = -2*rho[t]*rho[m]
end
if abs(m-t)<nt
aux += rho[abs(m-t)+1]
end
if m+t-1<=nt
aux += rho[m+t-1]
for t in 1:nt
is = max(1, t-iw-2) + 1
ie = is + iw-1
for m in is:ie
aux=0.0
if m <= nt
aux = -2*rho[t]*rho[m]
end
if abs(m-t)<nt
aux += rho[abs(m-t)+1]
end
if m+t-1<=nt
aux += rho[m+t-1]
end
dr[t] +=aux^2
end
dr[t] +=aux^2
end
dr[t] = sqrt(dr[t]/nd)
dr[t] = sqrt(dr[t]/nd)
end
return dr
end
end
"""
"""
bias(a1::ADerrors.uwreal [,a2::ADerrors.uwreal], id::Union{Int64,String} [, ws::Aderrors.wspace]; iw::Int64=0)
it return the bias associated to `id` that ADErrors includes in the estimate of the autocorrelation function
......@@ -137,97 +137,97 @@ function calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerro
if `iw = 0` the integration window is determined by the automatic windowing procedure with S_τ = 4.0
"""
function bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerrors.wspace;iw::Int64=0)
"""
function bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerrors.wspace;iw::Int64=0)
idx = ws.map_ids[id]
nd = ws.fluc[idx].nd
if nd ==1
return 0.0
return 0.0
end
gamma = calc_gamma(a1,a2,id,ws)
return gamma[1] + 2.0*sum(gamma[2:iw])
end
end
bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String,ws::ADerrors.wspace; iw::Int64=0)= bias(a1,a2,ws.str2id[id],ws,iw=iw)
bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64; iw::Int64=0) = bias(a1,a2,id,wsg,iw=iw)
bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String; iw::Int64=0)= bias(a1,a2,ADerrors.wsg.str2id[id],wsg,iw=iw)
bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String,ws::ADerrors.wspace; iw::Int64=0)= bias(a1,a2,ws.str2id[id],ws,iw=iw)
bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64; iw::Int64=0) = bias(a1,a2,id,wsg,iw=iw)
bias(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::String; iw::Int64=0)= bias(a1,a2,ADerrors.wsg.str2id[id],wsg,iw=iw)
function bias(a::ADerrors.uwreal,id::Int64,ws::ADerrors.wspace;iw::Int64=0)
function bias(a::ADerrors.uwreal,id::Int64,ws::ADerrors.wspace;iw::Int64=0)
idx = ws.map_ids[id]
nd = ws.fluc[idx].nd
if nd ==1
return 0.0
return 0.0
end
gamma = calc_gamma(a,id,ws)
iw = iw ==0 ? ADerrors.wopt_ulli(nd,ADerrors.DEFAULT_STAU,gamma) : iw
return gamma[1] + 2.0*sum(gamma[2:iw])
end
end
bias(a::ADerrors.uwreal,id::String,ws::ADerrors.wspace; iw::Int64=0) = bias(a,ws.str2id[id],ws,iw=iw)
bias(a::ADerrors.uwreal,id::Int64; iw::Int64=0) = bias(a,id,ADerrors.wsg,iw=iw)
bias(a::ADerrors.uwreal,id::String; iw::Int64=0) = bias(a,ADerrors.wsg.str2id[id],ADerrors.wsg,iw=iw)
bias(a::ADerrors.uwreal,id::String,ws::ADerrors.wspace; iw::Int64=0) = bias(a,ws.str2id[id],ws,iw=iw)
bias(a::ADerrors.uwreal,id::Int64; iw::Int64=0) = bias(a,id,ADerrors.wsg,iw=iw)
bias(a::ADerrors.uwreal,id::String; iw::Int64=0) = bias(a,ADerrors.wsg.str2id[id],ADerrors.wsg,iw=iw)
function bias(vobs::Vector{ADerrors.uwreal}, id::Int64, ws::ADerrors.wspace;)
function bias(vobs::Vector{ADerrors.uwreal}, id::Int64, ws::ADerrors.wspace;)
b = zeros(length(vobs),length(vobs))
iw = maximum([ADerrors.window(obs,id) for obs in vobs])
for i in 1:length(vobs)-1, j in i+1:length(vobs)
b[i,j] = bias(vobs[i],vobs[j],id,ws,iw=iw)
b[j,i] = b[i,j]
b[i,j] = bias(vobs[i],vobs[j],id,ws,iw=iw)
b[j,i] = b[i,j]
end
for i in eachindex(vobs)
b[i,i] = bias(vobs[i],id,ws,iw = 0)
b[i,i] = bias(vobs[i],id,ws,iw = 0)
end
return b
end
end
bias(vobs::Vector{ADerrors.uwreal}, id::String, ws::ADerrors.wspace) = bias(vobs,ws.str2id[id],ws)
bias(vobs::Vector{ADerrors.uwreal}, id::Int64) = bias(vobs,id,ADerrors.wsg)
bias(vobs::Vector{ADerrors.uwreal}, id::String) = bias(vobs,ADerrors.wsg.str2id[id],ADerrors.wsg)
bias(vobs::Vector{ADerrors.uwreal}, id::String, ws::ADerrors.wspace) = bias(vobs,ws.str2id[id],ws)
bias(vobs::Vector{ADerrors.uwreal}, id::Int64) = bias(vobs,id,ADerrors.wsg)
bias(vobs::Vector{ADerrors.uwreal}, id::String) = bias(vobs,ADerrors.wsg.str2id[id],ADerrors.wsg)
@doc raw"""
@doc raw"""
cov_pe(vobs::Vectir{ADerrors.uwreal}[,ws:ADerrors.wspace])
Computes the unbiased covariance matrix of `vobs` à la `pyerrors`.
It first computes the naive correlation matrix and the it multiply it with the errors
of the `vobs`. This method is exact if
of the `vobs`. This method is exact if
tau_int^ij = sqrt(tau_int^i tau_int^j)
It assumes that `ADerrors.uwerr` has been run on `vobs`
This method is stable in the sense that the covarince of two `ADerrors.uwreal` depends only
This method is stable in the sense that the covarince of two `ADerrors.uwreal` depends only
on them and not on the others `ADerrors.uwreal` that are in `vobs`
See also [`cov_min`](@ref) [`cov_ad`](@ref)
"""
function cov_pe(vobs::AbstractVector{ADerrors.uwreal}, ws::ADerrors.wspace)
function cov_pe(vobs::Vector{ADerrors.uwreal}, ws::ADerrors.wspace)
ids = ADerrors.unique_ids_multi(vobs, ws)
nid = length(ids)
nobs = length(vobs)
cov = zeros(Float64, nobs,nobs)
# naive covarince matrix
# naive covarince matrix
for j in 1:nid
idx = ws.map_ids[ids[j]]
nd_eff = div(ws.fluc[idx].nd, ws.fluc[idx].ibn)
for n1 in 1:nobs, n2 in n1:nobs
gamma = calc_gamma(vobs[n1],vobs[n2],ids[j],ws)
cov[n1,n2] +=gamma[1]/nd_eff
end
end
idx = ws.map_ids[ids[j]]
nd_eff = div(ws.fluc[idx].nd, ws.fluc[idx].ibn)
for n1 in 1:nobs, n2 in n1:nobs
gamma = calc_gamma(vobs[n1],vobs[n2],ids[j],ws)
cov[n1,n2] +=gamma[1]/nd_eff
end
end
for n1 in 1:nobs-1, n2 in n1+1:nobs ## symmetrizing
cov[n2,n1] = cov[n1,n2]
cov[n2,n1] = cov[n1,n2]
end
corr = [cov[i,j]/sqrt(cov[i,i]*cov[j,j]) for i in 1:nobs, j in 1:nobs]
cov = [corr[i,j]*vobs[i].err*vobs[j].err for i in 1:nobs, j in 1:nobs]
return cov
end
end
cov_pe(vobs::AbstractVector{ADerrors.uwreal}) = cov_pe(vobs,ADerrors.wsg)
cov_pe(vobs::Vector{ADerrors.uwreal}) = cov_pe(vobs,ADerrors.wsg)
@doc raw"""
@doc raw"""
cov_min(vobs::Vector{ADerrors.uwreal} [ws::ADerrors.wspace])
It computes the unbiased covariance matrix of `vobs` using a mixed method.
......@@ -237,55 +237,55 @@ function calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerro
matrix "à la pyerrors, that is, it multipy the covarince matrix with the errors of `vobs`
It assumes that `ADerrors.uwerr` has been run on `vobs`
This method is not stable, in the sense that the covarince between two `ADerrors.uwreal`
This method is not stable, in the sense that the covarince between two `ADerrors.uwreal`
depend also on the other `ADerrors.uwreal` included in `vobs`
See also [`cov_pe`](@ref) [`cov_ad`](@ref)
"""
function cov_min(vobs::AbstractVector{ADerrors.uwreal},ws::ADerrors.wspace)
function cov_min(vobs::Vector{ADerrors.uwreal},ws::ADerrors.wspace)
ids = ADerrors.unique_ids_multi(vobs, ws)
nid = length(ids)
iw = fill(typemax(Int64),nid)
iw = fill(typemax(Int64),nid)
for obs in vobs, j in eachindex(obs.ids), i in eachindex(ids)
if (obs.ids[j] == ids[i])
iw[i] = min(iw[i], obs.cfd[j].iw)
end
if (obs.ids[j] == ids[i])
iw[i] = min(iw[i], obs.cfd[j].iw)
end
end
nobs = length(vobs)
cov = zeros(Float64, nobs,nobs)
for j in 1:nid
idx = ws.map_ids[ids[j]]
nd_eff = div(ws.fluc[idx].nd, ws.fluc[idx].ibn)
for n in eachindex(vobs)
gamma = calc_gamma(vobs[n],ids[j],ws)
cov[n,n] += length(gamma)==1 ? gamma[1]/nd_eff : (gamma[1] + 2.0*sum(gamma[2:iw[j]]))/nd_eff
end
for n1 in 1:nobs-1, n2 in n1+1:nobs
gamma12 = calc_gamma(vobs[n1],vobs[n2],ids[j],ws)
if length(gamma12) == 1
cov[n1,n2] +=gamma12[1]/nd_eff
continue;
idx = ws.map_ids[ids[j]]
nd_eff = div(ws.fluc[idx].nd, ws.fluc[idx].ibn)
for n in eachindex(vobs)
gamma = calc_gamma(vobs[n],ids[j],ws)
cov[n,n] += length(gamma)==1 ? gamma[1]/nd_eff : (gamma[1] + 2.0*sum(gamma[2:iw[j]]))/nd_eff
end
for n1 in 1:nobs-1, n2 in n1+1:nobs
gamma12 = calc_gamma(vobs[n1],vobs[n2],ids[j],ws)
if length(gamma12) == 1
cov[n1,n2] +=gamma12[1]/nd_eff
continue;
end
gamma21 = calc_gamma(vobs[n2],vobs[n1],ids[j],ws)
cov[n1,n2] += (gamma12[1] + sum(gamma12[2:iw[j]]) + sum(gamma21[2:iw[j]]))/nd_eff
end
gamma21 = calc_gamma(vobs[n2],vobs[n1],ids[j],ws)
cov[n1,n2] += (gamma12[1] + sum(gamma12[2:iw[j]]) + sum(gamma21[2:iw[j]]))/nd_eff
end
end
for n1 in 1:nobs-1, n2 in n1+1:nobs
cov[n2,n1] = cov[n1,n2]
cov[n2,n1] = cov[n1,n2]
end
corr = [cov[i,j]/sqrt(cov[i,i]*cov[j,j]) for i in 1:nobs, j in 1:nobs]
cov = [corr[i,j]*vobs[i].err*vobs[j].err for i in 1:nobs, j in 1:nobs]
return cov
end
end
cov_min(vobs::AbstractVector{ADerrors.uwreal}) = cov_min(vobs,ADerrors.wsg)
cov_min(vobs::Vector{ADerrors.uwreal}) = cov_min(vobs,ADerrors.wsg)
@doc raw"""
@doc raw"""
cov_AD(vobs::Vector{ADerrors.uwrela} [ws::ADerrors.wspace]; biased::Bool=true, info::Bool=false)
It computes the unbiased covariance matrix of `vobs` "à la ADerrors".
......@@ -294,116 +294,116 @@ function calc_gamma(a1::ADerrors.uwreal,a2::ADerrors.uwreal,id::Int64,ws::ADerro
and computes the covarince matrix integrating the autocorrelation function up to `IW`
It assumes that `ADerrors.uwerr` has been run on `vobs`
This method is not stable, in the sense that the covarince between two `ADerrors.uwreal`
This method is not stable, in the sense that the covarince between two `ADerrors.uwreal`
depend also on the other `ADerrors.uwreal` included in `vobs`
# KEYWORD ARGUMENTS
- `bias::Bool`: if `true` (default), it computes the covariance matrix including the bias as `ADerrors`
- `info::Bool`: if `true` return the integration window used to compute the covariance matrix as a Dict{Int64,Int64}(id=> iw)
See also [`cov_min`](@ref) [`cov_pe`](@ref)
"""
function cov_AD(vobs::AbstractVector{ADerrors.uwreal},ws::ADerrors.wspace; biased::Bool=true, info::Bool=false)
function cov_AD(vobs::Vector{ADerrors.uwreal},ws::ADerrors.wspace; biased::Bool=true, info::Bool=false)
ids = ADerrors.unique_ids_multi(vobs, ws)
nid = length(ids)
iw = zeros(Int64,nid)
iw = zeros(Int64,nid)
for obs in vobs, j in eachindex(obs.ids), i in eachindex(ids)
if (obs.ids[j] == ids[i])
iw[i] = max(iw[i], obs.cfd[j].iw)
end
if (obs.ids[j] == ids[i])
iw[i] = max(iw[i], obs.cfd[j].iw)
end
end
nobs = length(vobs)
cov = zeros(Float64, nobs,nobs)
for j in 1:nid
idx = ws.map_ids[ids[j]]
nd_eff = div(ws.fluc[idx].nd, ws.fluc[idx].ibn)
for n in eachindex(vobs)
gamma = calc_gamma(vobs[n],ids[j],ws)
if biased
gamma .+= bias(vobs[n],ids[j],ws)/nd_eff
end
cov[n,n] += length(gamma)==1 ? gamma[1]/nd_eff : (gamma[1] + 2*sum(gamma[2:iw[j]]))/nd_eff
end
for n1 in 1:nobs-1, n2 in n1+1:nobs
gamma12 = calc_gamma(vobs[n1],vobs[n2],ids[j],ws)
if length(gamma12) == 1
cov[n1,n2] += gamma12[1]/nd_eff
continue;
idx = ws.map_ids[ids[j]]
nd_eff = div(ws.fluc[idx].nd, ws.fluc[idx].ibn)
for n in eachindex(vobs)
gamma = calc_gamma(vobs[n],ids[j],ws)
if biased
gamma .+= bias(vobs[n],ids[j],ws)/nd_eff
end
cov[n,n] += length(gamma)==1 ? gamma[1]/nd_eff : (gamma[1] + 2*sum(gamma[2:iw[j]]))/nd_eff
end
gamma21 = calc_gamma(vobs[n2],vobs[n1],ids[j],ws)
if biased
cov[n1,n2] += (gamma12[1] + sum(gamma12[2:iw[j]]) + sum(gamma21[2:iw[j]]))/nd_eff*(1+(2*iw[j]-1)/nd_eff)
else
cov[n1,n2] += (gamma12[1] + sum(gamma12[2:iw[j]]) + sum(gamma21[2:iw[j]]))/nd_eff
for n1 in 1:nobs-1, n2 in n1+1:nobs
gamma12 = calc_gamma(vobs[n1],vobs[n2],ids[j],ws)
if length(gamma12) == 1
cov[n1,n2] += gamma12[1]/nd_eff
continue;
end
gamma21 = calc_gamma(vobs[n2],vobs[n1],ids[j],ws)
if biased
cov[n1,n2] += (gamma12[1] + sum(gamma12[2:iw[j]]) + sum(gamma21[2:iw[j]]))/nd_eff*(1+(2*iw[j]-1)/nd_eff)
else
cov[n1,n2] += (gamma12[1] + sum(gamma12[2:iw[j]]) + sum(gamma21[2:iw[j]]))/nd_eff
end
end
end
end
for n1 in 1:nobs-1, n2 in n1+1:nobs
cov[n2,n1] = cov[n1,n2]
cov[n2,n1] = cov[n1,n2]
end
return info ? (cov, Dict(ids.=>iw)) : cov
end
end
cov_AD(vobs::AbstractVector{ADerrors.uwreal};biased::Bool=false) = cov_AD(vobs,ADerrors.wsg,biased=biased)
cov_AD(vobs::Vector{ADerrors.uwreal};biased::Bool=false) = cov_AD(vobs,ADerrors.wsg,biased=biased)
"""
"""
chiexp(hess, C, W)
It computes the expected chisquare given the hessian of the chisquare function computed at its minimum,
the covariance matrix of data and the weigth function.
It computes the expected chisquare given the hessian of the chisquare function computed at its minimum,
the covariance matrix of data and the weigth function.
"""
function chiexp(hess,C,W)
function chiexp(hess,C,W)
N, = size(hess)
ndata, = size(C)
npar = N-ndata;
L = LinearAlgebra.cholesky(W).L
Linv = LinearAlgebra.pinv(L)
_hess = view(hess, 1:npar,npar+1:npar+ndata)
S = [sum(Linv[alpha,gamma]*_hess[i,gamma] for gamma in 1:ndata) for i in 1:npar, alpha in 1:ndata]
H = [sum(S[n,alpha]*S[m,alpha] for alpha in 1:ndata) for n in 1:npar, m in 1:npar]
Hinv = LinearAlgebra.pinv(H)
P = [ sum(_hess[i,alpha]*Hinv[i,j]*_hess[j,beta] for i in 1:npar, j in 1:npar) for alpha in 1:ndata, beta in 1:ndata]
chiexp = let aux = C*(W-P)
sum(aux[alpha,alpha] for alpha in 1:ndata )
sum(aux[alpha,alpha] for alpha in 1:ndata )
end
return chiexp
end
"""
return chiexp
end
"""
chiexp(chisq::Function,par,d,C,W)
It computes the expected chisquare given the chisquare function, the fit parameters and the data.
This is preferred to ADerrors.chiexp for correlated fits since it accepts the covariance matrix as an input,
It computes the expected chisquare given the chisquare function, the fit parameters and the data.
This is preferred to ADerrors.chiexp for correlated fits since it accepts the covariance matrix as an input,
that usually has been already computed before hand.
"""
function chiexp(chisq::Function,par::AbstractVector{Float64},d::AbstractVector{Float64},C::AbstractMatrix{Float64},W::AbstractMatrix{Float64})
function chiexp(chisq::Function,par::Vector{Float64},d::Vector{Float64},C::AbstractMatrix{Float64},W::AbstractMatrix{Float64})
ndata = length(d);
if size(W) != (ndata,ndata)
throw(DimensionMismatch("[chiexp] d and W must have same dimension"))
throw(DimensionMismatch("[chiexp] d and W must have same dimension"))
end
npar = length(par);
hess = zeros(npar+ndata,npar+ndata);
ccsq(x) = chisq(view(x,1:npar),view(x,npar+1:npar+ndata))
ForwardDiff.hessian!(hess,ccsq,[par; d])
return chiexp(hess,C,W)
end
end
function chiexp(chisq::Function,par::AbstractVector,d::AbstractVector{Float64},C::AbstractMatrix{Float64},W::AbstractVector{Float64})
function chiexp(chisq::Function,par::Vector,d::Vector{Float64},C::AbstractMatrix{Float64},W::Vector{Float64})
ndata = length(d);
if length(W) != ndata
throw(DimensionMismatch("[chiexp] d and W must have same dimension"))
throw(DimensionMismatch("[chiexp] d and W must have same dimension"))
end
npar = length(par);
hess = zeros(npar+ndata,npar+ndata);
ccsq(x) = chisq(view(x,1:npar),view(x,npar+1:npar+ndata))
ForwardDiff.hessian!(hess,ccsq,[par; d])
return chiexp(hess,C,LinearAlgebra.Diagonal(W))
end
end
chiexp(chisq::Function,par::AbstractVector{ADerrors.uwreal},d::AbstractVector{ADerrors.uwreal},C,W) = chiexp(chisq,ADerrors.value.(par),ADerrors.value.(d),C,W)
chiexp(chisq::Function,par::Vector{ADerrors.uwreal},d::Vector{ADerrors.uwreal},C,W) = chiexp(chisq,ADerrors.value.(par),ADerrors.value.(d),C,W)
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment