Returns an `uwreal` data type. Depending on the first argument, the `uwreal` stores the following information:
- A single `Float64`. In this case the variable acts in exactly the same way as a real number. This is understood as a quantity with zero error.
- A 2 element Vector of `Float64` `[value, error]`. In this case the data is understood as `value +/- error`.
- A Vector of `Float64` of length larger than 2. In this case the data is understood as consecutive measurements of an observable in a Monte Carlo (MC) simulation.
- Input is a single `Float64`. In this case the variable acts in exactly the same way as a real number. This is understood as a quantity with zero error.
- Input is a 2 element Vector of `Float64` `[value, error]`. In this case the data is understood as `value +/- error`.
- Input is a Vector of `Float64` of length larger than 4. In this case the data is understood as consecutive measurements of an observable in a Monte Carlo (MC) simulation.
In the last two cases, an ensemble `ID` is required as input. Data with the same `ID` are considered as correlated (i.e. fully correlated for the case of a `value +/- error` observables measured on the same sample for the case of data from a MC simulation). For example:
In some situations an observable is not measured in every configuration. In this case two additional arguments aree needed to define the observable
- `idm`. Type `Vector{Int64}`. `idm[n]` labels the configuration where data[n] is measured.
In some situations an observable is not measured in every configuration. In this case two additional arguments are needed to define the observable
- `idm`. Type `Vector{Int64}`. `idm[n]` labels the configuration where `data[n]`is measured.
- `nms`. Type `Int64`. The total number of measurements in the ensemble
```@example
using ADerrors # hide
...
...
@@ -694,17 +703,17 @@ b = sin(2.0*a)
uwerr(b)
println("LookhowIpropagateerrors:", b)
c = 1.0 + b - 2.0*sin(a)*cos(a)
uwerr(b)
uwerr(c)
println("LookhowgoodIamatthis(zeroerror!):", c)
```
### Optimal window
Error in data coming from a Monte Carlo ensemble is performed by summing the autocorrelation function ``\Gamma_{\rm ID}(t)`` of the data for each ensemble ID. In practice this sum is truncated up to a window ``W_{\rm ID}``.
Error in data coming from a Monte Carlo ensemble is determined by summing the autocorrelation function ``\Gamma_{\rm ID}(t)`` of the data for each ensemble ID. In practice this sum is truncated up to a window ``W_{\rm ID}``.
By default, the summation window is determined as [proposed by U. Wolff](https://inspirehep.net/literature/621085) with a parameter ``S_{\tau} = 4``, but other methods are avilable via the optional argument `wpm`.
For each ensemble `ID` one has to pass a vector of `Float64` of length 4. The first three components of the vector specify the criteria to determine the summation window:
For each ensemble `ID` one can pass a vector of `Float64` of length 4. The first three components of the vector specify the criteria to determine the summation window:
- `vp[1]`: The autocorrelation function is summed up to `t = round(vp[1])`.
- `vp[2]`: The sumation window is determined [using U. Wolff poposal](https://inspirehep.net/literature/621085) with ``S_{\tau} = {\rm vp[2]}``.
- `vp[3]`: The autocorrelation function ``\Gamma(t)`` is summed up a point where its error ``\delta\Gamma(t)`` is a factor `vp[3]` times larger than the signal.
...
...
@@ -716,44 +725,80 @@ Note that:
- One, and only one, of the components `vp[1:3]` has to be positive. This chooses your criteria to determine the summation window.
```@example
using ADerrors # hide
a = uwreal(rand(2000), 1233)
# Generate some correlated data
eta = randn(1000)
x = Vector{Float64}(undef, 1000)
x[1] = 0.0
for i in 2:1000
x[i] = x[i-1] + eta[i]
if abs(x[i]) > 1.0
x[i] = x[i-1]
end
end
# Load the data in a uwreal
a = uwreal(x.^2, 1233)
wpm = Dict{Int64,Vector{Float64}}()
# Use default analysis (stau = 4.0)
uwerr(a)
println(a)
println("default:", a, "(tauint=", taui(a, 1233), ")")
# This will still do default analysis because
# a does not depend on emsemble 300
wpm[300] = [-1.0, 8.0, -1.0, 145.0]
uwerr(a, wpm)
println(a)
println("default:", a, "(tauint=", taui(a, 1233), ")")
# Fix the summation window to 1 (i.e. uncorrelated data)
wpm[1233] = [1.0, -1.0, -1.0, -1.0]
uwerr(a, wpm)
println(a)
println("uncorrelated:", a, "(tauint=", taui(a, 1233), ")")
# Use stau = 1.5
wpm[1233] = [-1.0, 1.5, -1.0, -1.0]
uwerr(a, wpm)
println(a)
println("stau=1.5:", a, "(tauint=", taui(a, 1233), ")")
# Use fixed window 15 and add tail with texp = 100.0
wpm[1233] = [15.0, -1.0, -1.0, 100.0]
uwerr(a, wpm)
println(a)
println("Fixedwindow15,texp=100:", a, "(tauint=", taui(a, 1233), ")")
# Sum up to the point that the signal in Gamma is
# 1.5 times the error and add a tail with texp = 30.0
wpm[1233] = [-1.0, -1.0, 1.5, 30.0]
uwerr(a, wpm)
println(a)
println("signal/noise=1.5,texp=30:", a, "(tauint=", taui(a, 1233), ")")
An optional parameter `wpm` can be used to choose the summation window for the relevant autocorrelation functions. The situation is completely analogous to the case of error analysis of single variables.
@@ -11,53 +11,214 @@ function find_mcid(a::uwreal, mcid::Int64)
end
end
return0
returnnothing
end
err(a::uwreal)=a.err
"""
err(a::uwreal)
Returns the error of the `uwreal` variable `a`. It is assumed that `uwerr` has been run on the variable so that an error is available. Otherwise an error message is printed.
```@example
using ADerrors # hide
a = uwreal([1.2, 0.2], 12) # a = 1.2 +/- 0.2
uwerr(a)
println("ahaserror:", err(a))
```
"""
function err(a::uwreal)
if(length(a.cfd)==0)
error("No error available... maybe run uwerr")
end
returna.err
end
"""
value(a::uwreal)
Returns the (mean) value of the `uwreal` variable `a`
```@example
using ADerrors # hide
a = uwreal([1.2, 0.2], 12) # a = 1.2 +/- 0.2
uwerr(a)
println("ahascentralvalue:", value(a))
```
"""
value(a::uwreal)=a.mean
derror(a::uwreal)=a.derr
"""
derror(a::uwreal)
Returns an estimate of teh error of the error of the `uwreal` variable `a`. It is assumed that `uwerr` has been run on the variable so that an error is available. Otherwise an error message is printed.
```@example
using ADerrors # hide
a = uwreal([1.2, 0.2], 12) # a = 1.2 +/- 0.2
uwerr(a)
println("ahaserroroftheerror:", derror(a))
```
"""
function derror(a::uwreal)
if(length(a.cfd)==0)
error("No error available... maybe run uwerr")
end
returna.derr
end
"""
taui(a::uwreal, id::Int64)
Returns the value of tauint for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
eta = randn(1000)
x = Vector{Float64}(undef, 1000)
x[1] = 0.0
for i in 2:1000
x[i] = x[i-1] + eta[i]
if abs(x[i]) > 1.0
x[i] = x[i-1]
end
end
a = uwreal(x.^2, 666)
uwerr(a)
println("Erroranalysisresult:", a, "(tauint=", taui(a, 666), ")")
```
"""
function taui(a::uwreal,mcid::Int64)
idx=find_mcid(a,mcid)
if(idx==0)
return0.5
if(idx==nothing)
error("No error available... maybe run uwerr")
else
returna.cfd[idx].taui
end
end
"""
dtaui(a::uwreal, id::Int64)
Returns an estimate on the error of tauint for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
Returns the summation window for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
eta = randn(1000)
x = Vector{Float64}(undef, 1000)
x[1] = 0.0
for i in 2:1000
x[i] = x[i-1] + eta[i]
if abs(x[i]) > 1.0
x[i] = x[i-1]
end
end
a = uwreal(x.^2, 666)
uwerr(a)
println("Erroranalysisresult:", a,
"(window=", window(a, 666), ")")
```
"""
function window(a::uwreal,mcid::Int64)
idx=find_mcid(a,mcid)
if(idx==0)
return0
if(idx==nothing)
error("No error available... maybe run uwerr")
else
returna.cfd[idx].iw
end
end
"""
rho(a::uwreal, id::Int64)
Returns the normalized autocorrelation function of `a` for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
eta = randn(1000)
x = Vector{Float64}(undef, 1000)
x[1] = 0.0
for i in 2:1000
x[i] = x[i-1] + eta[i]
if abs(x[i]) > 1.0
x[i] = x[i-1]
end
end
a = uwreal(x.^2, 666)
uwerr(a)
v = rho(a, 666)
for i in 1:length(v)
println(i, "", v[i])
end
```
"""
function rho(a::uwreal,mcid::Int64)
idx=find_mcid(a,mcid)
if(idx==0)
return[1.0]
if(idx==nothing)
error("No error available... maybe run uwerr")
else
returna.cfd[idx].gamm./a.cfd[idx].gamm[1]
end
end
"""
rho(a::uwreal, id::Int64)
Returns an estimate of the error on the normalized autocorrelation function of `a` for the ensemble `id`. It is assumed that `uwerr` has been run on the variable and that `id` contributes to the observable `a`. Otherwise an error message is printed.
```@example
using ADerrors # hide
# Generate some correlated data
eta = randn(1000)
x = Vector{Float64}(undef, 1000)
x[1] = 0.0
for i in 2:1000
x[i] = x[i-1] + eta[i]
if abs(x[i]) > 1.0
x[i] = x[i-1]
end
end
a = uwreal(x.^2, 666)
uwerr(a)
v = rho(a, 666)
dv = drho(a, 666)
for i in 1:length(v)
println(i, "", v[i], "+/-", dv[i])
end
```
"""
function drho(a::uwreal,mcid::Int64)
idx=find_mcid(a,mcid)
if(idx==0)
return[0.0]
if(idx==nothing)
error("No error available... maybe run uwerr")
else
returna.cfd[idx].drho
end
...
...
@@ -219,11 +380,11 @@ function details(a::uwreal, ws::wspace, io::IO=stdout, names::Dict{Int64, String
Returns a vector of `uwreal` such that their mean values are `avgs[:]` and their covariance is `Mcov[:,:]`. In order to construct these observables `n=length(avgs)` ensemble ID are used. These have to be specified in the vector `ids[:]`.
```@example
using ADerrors # hide
# Put some average values and covariance
avg = [16.26, 0.12, -0.0038]
Mcov = [0.478071 -0.176116 0.0135305
-0.176116 0.0696489 -0.00554431
0.0135305 -0.00554431 0.000454180]
# Produce observables with ensemble ID
# [1, 2001, 32]. Do error analysis
p = cobs(avg, Mcov, [1, 2001, 32])
uwerr.(p)
# Check central values are ok
avg2 = value.(p)
println("Betterbezero:", sum((avg.-avg2).^2))
# Check that the covariance is ok
Mcov2 = cov(p)
println("Betterbezero:", sum((Mcov.-Mcov2).^2))
```
"""
function cobs(avgs::Vector{Float64},cov::Array{Float64,2},ids::Vector{Int64})
where the function ``f_i(p)`` is an arbitrary function of the fit parameters. In simple words, the function is assumed to be quadratic in the data.
- `xp`: A vector of `Float64`. The value of the fit parameters at the minima.
- `data`: A vector of `uwreal`. The data whose fluctuations enter in the evaluation of the `chisq`.
- `W`: A matrix. The weights that enter in the evaluation of the `chisq` function. If a vector is passed, the matrix is assumed to be diagonal (i.e. **uncorrelated** fit). If no weights are passed, the routines assumes that `W` is diagonal with entries given by the inverse errors squared of the data (i.w. the `chisq` is weighted with the errors of the data).
