API

uwreal(x::Float64)
uwreal([value::Float64, error::Float64], mcid)
uwreal(data::Vector{Float64}, mcid[, replica::Vector{Int64}])
uwreal(data::Vector{Float64}, mcid[, replica::Vector{Int64}], idm::Vector{Int64}, nms::Int64)

Returns an uwreal data type. Depending on the first argument, the uwreal stores the following information:

• 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). The preferred way to input the ensemble tag is via a String that uniquely identifies the ensemble, but an integer is also supported for legacy reasons. For example:

using ADerrors # hide
a = uwreal([1.2, 0.2], "Simple var with error")   # a = 1.2 +/- 0.2
b = uwreal([1.2, 0.2], "Simple var with error")   # b = 1.2 +/- 0.2
c = uwreal([1.2, 0.2], "Another var with error")  # c = 1.2 +/- 0.2

d = a-b
uwerr(d)
println("d has zero error because a and b are correlated", d)

e = a-c
uwerr(e)
println("e has non zero error because a and c are independent", e)

Replica

data can contain measurements in several replica (i.e. independent simulations with the same physical and algorithmic parameters). How many measurements correspond to each replica are specified by an optional replica vector argument replica. Note that length(replica) is the number of replica (i.e. independent simulations), and that sum(replica) must match length(data)

using ADerrors # hide
# 1000 measurements in three replica of lengths
# 500, 100 and 400
a = uwreal(rand(1000), "Ensemble with three replica", [500, 100, 400]) 

Gaps in the measurements

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

using ADerrors # hide
# Observable measured on the odd configurations 
# 1, 3, 5, ..., 999 on an emsemble of length 1000
a = uwreal(rand(500), "Observable with gaps", collect(1:2:999), 1000)

# Observable measured on the first 900 configurations 
# on the same emsemble
b = uwreal(rand(900), "Observable with gaps", collect(1:900), 1000)

Note that in this case, if the ensemble has different replica, sum(replica) must match nms

using ADerrors # hide
# Observable measured on the even configurations 
# 2, 4, 6, ..., 200 on an emsemble of length 200
# with two replica of lengths 75, 125
a = uwreal(data_a[1:500], "Observable with gaps in an ensemble with replica", [75, 125], collect(2:2:200), 200)

cobs(avgs::Vector{Float64}, Mcov::Array{Float64, 2}, str::String)
cobs(avgs::Vector{Float64}, Mcov::Array{Float64, 2}, ids::Vector{Int64})

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 are generated either from the string str or have to be specified in the vector ids[:].

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, "GF beta function parameters")
uwerr.(p)

# Check central values are ok
avg2 = value.(p)
println("Better be zero: ", sum((avg.-avg2).^2))

# Check that the covariance is ok
Mcov2 = cov(p)
println("Better be zero: ", sum((Mcov.-Mcov2).^2))

change_id(; from::String, to::String)

Changes the ensemble id from the string/integer from to the string to. Note that a call to this routine will affect the output of all currently defined uwreal types that depend on ensemble ID from. On the other hand, it will not affect data that is read after the call to change_id. This means that:

• All input should be performed before any calls to change_id.
• All output (i.e. calls to write_uwreal or details) is affected by previous calls to change_id

using ADerrors # hide
a = uwreal([1.2, 0.2], "Simple var with error")   # a = 1.2 +/- 0.2
b = uwreal([1.2, 0.2], 1233)  # c = 1.2 +/- 0.2

d = a-b
uwerr(d)
details(d)

change_id(from=1233, to="Error source from experiment")
change_id(from="Simple var with error", to="Error source from simulations")
details(d)

ensembles(a::uwreal)

Returns the list of ensembles in lexicographic order contributing to the error of a as a Vector{String}.

using ADerrors # hide
using BDIO
a = uwreal(rand(2000), "Ensemble A")
b = uwreal(rand(2034), "Ensemble B")
c = sin(a+b)
for i in ensembles(c)
    println("Ensemble "*str*" contributing")
end

write_uwreal(p::uwreal, fb, iu::Int)

Given a BDIO file handler fb, this writes the observable p in a BDIO record with user info iu.

using ADerrors # hide
using BDIO
a = uwreal(rand(2000), 12)

# Create a BDIO file and write observable with user info 8.
fb = BDIO_open("/tmp/foo.bdio", "w", "Test file")
write(a, fb, 8)
BDIO_close(fb)

# Open the file and move to first record
fb = BDIO_open("/tmp/foo.bdio", "r")
BDIO_seek!(fb)

# Read observable
b = read_uwreal(fb)

# Check
c = a - b
uwerr(c)
println("Better be zero: ", c)

read_uwreal(fb[, map_ids::Dict{Int64, Int64}])

Given a BDIO file handler fb, this routine returns the observable stored in the current record. Optionally, ensemble ID stored in the file can be changed at read time by passing a dictionary map_ids

using ADerrors # hide
using BDIO
a = uwreal(rand(2000), 12)

fb = BDIO_open("/tmp/foo.bdio", "w", "Test file")
write(a, fb, 8)
BDIO_close(fb)

# Open the file and move to first record
fb = BDIO_open("/tmp/foo.bdio", "r")
BDIO_seek!(fb)

# Read observable
changeID_12_by_120 = Dict(12 => 120)
b = read_uwreal(fb, changeID_12_by_120)

# Check
c = a - b
uwerr(c)
print("Better be zero: ")
details(c)

 uwerr(a::uwreal[, wpm::Dict{Int64,Vector{Float64}}])
 uwerr(a::uwreal[, wpm::Dict{String,Vector{Float64}}])

Performs error analysis on the observable a.

using ADerrors # hide
a = uwreal([1.3, 0.01], "Var with error") # 1.3 +/- 0.01
b = sin(2.0*a)
uwerr(b)
println("Look how I propagate errors:              ", b)
c = 1.0 + b - 2.0*sin(a)*cos(a)
uwerr(c)
println("Look how good I am at this (zero error!): ", c)

Optimal window

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 with a parameter $S_{\tau} = 4$, but other methods are avilable via the optional argument wpm.

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 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.

An additional fourth parameter vp[4], tells ADerrors to add a tail to the error with $\tau_{\rm exp} = {\rm vp[4]}$. See the reference for an explanation of the procedure.

Note that:

• Negative values of vp[1:4] are ignored.
• One, and only one, of the components vp[1:3] has to be positive. This chooses your criteria to determine the summation window.

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

# Load the data in a uwreal
a = uwreal(x.^2, "Random walk in [-1,1]")
wpm = Dict{String,Vector{Float64}}()

# Use default analysis (stau = 4.0)
uwerr(a)
println("default:                   ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")

# This will still do default analysis because 
# a does not depend on emsemble foo
wpm["Ensemble foo"] = [-1.0, 8.0, -1.0, 145.0]
uwerr(a, wpm)
println("default:                   ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")

# Fix the summation window to 1 (i.e. uncorrelated data)
wpm["Random walk in [-1,1]"] = [1.0, -1.0, -1.0, -1.0]
uwerr(a, wpm)
println("uncorrelated:              ",  a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")

# Use stau = 1.5
wpm["Random walk in [-1,1]"] = [-1.0, 1.5, -1.0, -1.0]
uwerr(a, wpm)
println("stau = 1.5:                ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")

# Use fixed window 15 and add tail with texp = 100.0
wpm["Random walk in [-1,1]"] = [15.0, -1.0, -1.0, 100.0]
uwerr(a, wpm)
println("Fixed window 15, texp=100: ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")

# Sum up to the point that the signal in Gamma is 
# 1.5 times the error and add a tail with texp = 10.0
wpm["Random walk in [-1,1]"] = [-1.0, -1.0, 1.5, 30.0]
uwerr(a, wpm)
println("signal/noise=1.5, texp=10: ", a, " (tauint = ", taui(a, "Random walk in [-1,1]"), ")")

cov(a::Vector{uwreal}[, wpm])

Determine the covariance matrix between the vector of observables a[:].

using ADerrors, LinearAlgebra # hide
a = uwreal([1.3, 0.01], "Var 1") # 1.3 +/- 0.01
b = uwreal([5.3, 0.23], "Var 2") # 5.3 +/- 0.23
uwerr(a)
uwerr(b)

x = [a+b, a-b]
mat = ADerrors.cov(x)
println("Covariance: ", mat[1,1], " ", mat[1,2])
println("            ", mat[2,1], " ", mat[2,2])
println("Check (should be zero): ",  mat[1,1] - mat[2,2])
println("Check (should be zero): ",  mat[1,2] - (err(a)^2-err(b)^2))

Case of Monte Carlo data

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.

 trcov(M::Array{Float64, 2}, a::Vector{uwreal}[, wmp])

Given a vector of uwreal, a[:] and a two dimensional symmtric positive definite array M, this routine computes ${\rm tr}(MC)$, where $C_{ij} = {\rm cov}(a[i], a[j])$.

using ADerrors, LinearAlgebra # hide
a = uwreal([1.3, 0.01], "Var with error 1") # 1.3 +/- 0.01
b = uwreal([5.3, 0.23], "Var with error 2") # 5.3 +/- 0.23
c = uwreal(rand(2000), "White noise ensemble")

x = [a+b+sin(c), a-b+cos(c), c-b/a]
M = [1.0 0.2 0.1
     0.2 2.0 0.3
     0.1 0.3 1.0]

mcov = cov(x)
d = tr(mcov * M)
println("Better be zero: ", d -trcov(M, x))

# Case of Monte Carlo data

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.

value(a::uwreal)

Returns the (mean) value of the uwreal variable a

using ADerrors # hide
a = uwreal([1.2, 0.2], 12)   # a = 1.2 +/- 0.2
uwerr(a)
println("a has central value: ", value(a))

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.

using ADerrors # hide
a = uwreal([1.2, 0.2], 12)   # a = 1.2 +/- 0.2
uwerr(a)
println("a has error: ", err(a))

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.

using ADerrors # hide
a = uwreal([1.2, 0.2], 12)   # a = 1.2 +/- 0.2
uwerr(a)
println("a has error of the error: ", derror(a))

taui(a::uwreal, mcid)

Returns the value of tauint for the ensemble mcid. It is assumed that uwerr has been run on the variable and that mcid contributes to the observable a. Otherwise an error message is printed. mcid can be either an Int64 (the proper ensemble ID), or a String (the ensemble tag).

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, "Some simple ensemble")
uwerr(a)
println("Error analysis result: ", a, " (tauint = ", taui(a, "Some simple ensemble"), ")")

dtaui(a::uwreal, mcid)

Returns an estimate on the error of tauint for the ensemble mcid. It is assumed that uwerr has been run on the variable and that mcid contributes to the observable a. Otherwise an error message is printed.

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, "Some simple ensemble")
uwerr(a)
println("Error analysis result: ", a, 
        " (tauint = ", taui(a, "Some simple ensemble"), " +/- ", dtaui(a, "Some simple ensemble"), ")")

rho(a::uwreal, mcid)

Returns the normalized autocorrelation function of a for the ensemble mcid. It is assumed that uwerr has been run on the variable and that mcid contributes to the observable a. Otherwise an error message is printed.

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, "Some simple ensemble")
uwerr(a)
v = rho(a, "Some simple ensemble")
for i in 1:length(v)
    println(i, " ", v[i])
end

drho(a::uwreal, mcid)

Returns an estimate of the error on the normalized autocorrelation function of a for the ensemble mcid. It is assumed that uwerr has been run on the variable and that mcid contributes to the observable a. Otherwise an error message is printed.

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, "Some simple ensemble")
uwerr(a)
v  =  rho(a, "Some simple ensemble")
dv = drho(a, "Some simple ensemble")
for i in 1:length(v)
    println(i, " ", v[i], " +/- ", dv[i])
end

mchist(a::uwreal, mcid)

Returns the fluctuations of the uwreal variable a on ensemble mcid. It is assumed that uwerr has been run on the variable and that mcid contributes to the observable a. Otherwise an error message is printed.

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, "Some simple ensemble")
uwerr(a)
v = mchist(a, "Some simple ensemble")
for i in 1:length(v)
    println(i, " ", v[i])
end

window(a::uwreal, mcid)

Returns the summation window for the ensemble mcid. It is assumed that uwerr has been run on the variable and that mcid contributes to the observable a. Otherwise an error message is printed.

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, "Some simple ensemble")
uwerr(a)
println("Error analysis result: ", a, 
        " (window = ", window(a, "Some simple ensemble"), ")")

details(a::uwreal; io::IO=stdout, names::Dict{Int64, String} = Dict{Int64, String}())
details(a::uwreal, ensemble::String, ws::wspace, io::IO=stdout)

Write out a detailed information on the error of a.

Arguments

If a the ergument ensemble is present, this routine writes information of the fluctuations of the observable in this ensemble. Optionally one can pass as a keyword argument (io) the IO stream to write to.

Example

using ADerrors # hide
a = uwreal(rand(2000),   "Ensemble A12", ["A12 REP 0", "A12 REP 1", "A12 REP 2"], [1000, 500, 500])
b = uwreal([1.2, 0.023], "Ensemble XYZ")
c = uwreal([5.2, 0.03],  "Ensemble RRR")
d = a + b - c
uwerr(d)

details(d)
for en in ensembles(d)
    details(d, en)
end

neid(a::uwreal)

Returns the number of different ensemble ID's contributing to a

using ADerrors # hide
a = uwreal([1.2, 0.2], 12)   # a = 1.2 +/- 0.2
b = uwreal([7.2, 0.5], 13)   # a = 7.2 +/- 0.5

c = a*b
println("Number od ID contributing to c: ", neid(a))

derivative(a::uwreal, b::uwreal)

Determines the derivate of the observable a with respect to observable b (i.e. `rac{{ m d} a}{{ m d} b}). Observablebmust depend on only a single ensemble (typicallyb` is just a primary observable).

using ADerrors # hide

# Put some data 
a = uwreal(1.2 .+ randn(2000), "Ensemble test A", ["R0", "Rep 1", "Last rep"], [1000, 550, 450])
b = uwreal(0.8 .+ randn(2034), "Ensemble test B", [1000, 584, 450])

d = sin(a+b)
println("Derivative of d w.r.t a:      ", derivative(d, a))
println("Derivative of d w.r.t log(a): ", derivative(d, log(a)))

root_error(fnew::Function, x0::Float64, data::Vector{uwreal})

Returns x such that fnew(x, data) = 0 (i.e. a root of the function). The error in data is propagated to the root position x.

using ADerrors # hide
# First define some arbitrary data
data = Vector{uwreal}(undef, 3)
data[1] = uwreal([1.0, 0.2],   "Var A")
data[2] = uwreal([1.2, 0.023], "Var B")
data[3] = uwreal(rand(1000),   "White noise ensemble")

# Now define a function
f(x, p) = x + p[1]*x + cos(p[2]*x+p[3])

# Find its root using x0=1.0 as initial
# guess of the position of the root
x = root_error(f, 1.0, data)
uwerr(x)
println("Root: ", x)

# Check
z = f(x, data)
uwerr(z)
print("Better be zero (with zero error): ")
details(z)

fit_error(chisq::Function, xp::Vector{Float64}, data::Vector{uwreal}[, wpm];
               W = Vector{Float64}(), chi_exp = true)

Given a $\chi^2(p, d)$, function of the fit parameters p[:] and the data d[:], this routine return the fit parameters as uwreal type and optionally, the expected value of $\chi^2(p, d)$.

Arguments

• chisq: Must be a function of two vectors (i.e. chisq(p::Vector, d::Vector)). To determine the fit parameters, the function can be arbitrary, but for the determination of the expected $\chi^2(p, d)$ is assumed to have the form

$\chi^2(p, d) = \sum_{ij} [d_i - f_i(p)]W_{ij}[d_j - f_j(p)]$

where the function $f_i(p)$ is an arbitrary function of the fit parameters. In simple words, the expected $\chi^2(p, d)$ is determined assuming that the function $\chi^2(p, d)$ is 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.
• wpm: Dict{Int64,Vector{Float64}} or Dict{String,Vector{Float64}}. The criteria to determine the summation window. See the documentation on uwerr function for more details.
• 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.e. the chisq is weighted with the errors of the data).
• chi_exp: Bool type. If false, do not compute the expected $\chi^2(p, d)$.

Example

using ADerrors, Distributions # hide

# Generate correlated samples with average 0.1
npt = 12
sig = zeros(npt, npt)
dx  = zeros(npt)
for i in 1:npt
    dx[i]    = 0.01*i
    sig[i,i] = dx[i]^2
    for j in i+1:npt
        sig[i,j] = 0.0001 - 0.000005*abs(i-j)
        sig[j,i] = 0.0001 - 0.000005*abs(i-j)
    end
end
dmv = MvNormal([0.1 for n in 1:npt], sig)
vs  = rand(dmv, 1)

# Create the uwreal data that we want to 
# fit to a constant
dt = cobs(vs[:,1], sig, "Data points")

# Define the chi^2
chisq(p, d) = sum( (d .- p[1]) .^ 2 ./ dx .^2 )

# The result of an uncorrelated fit to a 
# constant is the weighted average
xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]


# Propagate errors to the fit parameters and
# determine the expected chi^2
(fitp, csqexp) = fit_error(chisq, xp, dt)

uwerr.(fitp)
println(" *** FIT RESULTS ***")
print("Fit parameter:     ")
details.(fitp)
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", csqexp, "  (dof: ", npt-1, ")")

chiexp(chisq::Function,
            xp::Vector{Float64}, 
            data::Vector{uwreal};
            W = Vector{Float64}())

Given a $\chi^2(p, d)$, function of the fit parameters p[:] and the data d[:], compute the expected value of the $\chi^2(p, d)$.

Arguments

• chisq: Must be a function of two vectors (i.e. chisq(p::Vector, d::Vector)). The function is assumed to have the form

$\chi^2(p, d) = \sum_i [d_i - f_i(p)]W_{ij}[d_j - f_j(p)]$

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).

Example

using ADerrors, Distributions # hide

# Generate correlated samples with average 0.1
npt = 12
sig = zeros(npt, npt)
dx  = zeros(npt)
for i in 1:npt
    dx[i]    = 0.01*i
    sig[i,i] = dx[i]^2
    for j in i+1:npt
        sig[i,j] = 0.0001 - 0.000005*abs(i-j)
        sig[j,i] = 0.0001 - 0.000005*abs(i-j)
    end
end
dmv = MvNormal([0.1 for n in 1:npt], sig)
vs  = rand(dmv, 1)

# Create the uwreal data that we want to 
# fit to a constant
dt = cobs(vs[:,1], sig, [100+n for n in 1:npt])

# Define the chi^2
chisq(p, d) = sum( (d .- p[1]) .^ 2 ./ dx .^2 )

# The result of an uncorrelated fit to a 
# constant is the weighted average
xp = [sum(value.(dt) ./ dx)/sum(1.0 ./ dx)]

# Compare chi^2 and expected chi^2
println("chi^2 / chi_exp^2: ", chisq(xp, value.(dt)), " / ", chiexp(chisq, xp, dt))

int_error(fint::Function, a, b, p::Vector{uwreal})

Computes the integral

$\int_a^b {\rm fint}(x; p)\, {\rm d}x$

The errors in the parameters p[:] and (optionally) in the limits of the intgeral (a, b) are propagated to an error of the integral. The result is returned as an uwreal.

using ADerrors

# Average values and covariance from
# https://inspirehep.net/literature/1477411
# here we try to reproduce the result
# Scale ratio = 21.86(42) Eq.5.6
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]

p = cobs(avg, Mcov, "Beta function fit parameters")
g1s = uwreal([2.6723, 0.0064], 4)
g2s = 11.31

fint(x, p) = - (p[1] + p[2]*x^2 + p[3]*x^4)/x^3
g1 = sqrt(g1s)
g2 = sqrt(g2s)
sint = 2.0*exp(-int_error(fint, g1, g2, p))
uwerr(sint)
print("  From integral evaluation: ")
details(sint)