`ADerrors.jl` is a package for error propagation and analysis of Monte carlo data. At the core of the package is the `uwreal` data type, that is able to store variables with uncertainties
`ADerrors.jl` is a package for error propagation and analysis of Monte
carlo data. At the core of the package is the `uwreal` data type. It
contains variables with uncertainties. In very simple cases this is
just a central value and an error
```@repl gs
```@repl gs
using ADerrors
using ADerrors
a = uwreal([1.0, 0.1], 1) # 1.0 +/- 0.1
a = uwreal([1.0, 0.1], 1) # 1.0 +/- 0.1
```
```
It can also store MC data
But it also supports the case where the uncertainty comes from Monte
Carlo measurements
```@repl gs
```@repl gs
# Generate some correlated data
# Generate some MC correlated data
eta = randn(1000);
eta = randn(1000);
x = Vector{Float64}(undef, 1000);
x = Vector{Float64}(undef, 1000);
x[1] = 0.0;
x[1] = 0.0;
...
@@ -22,7 +26,7 @@ end
...
@@ -22,7 +26,7 @@ end
b = uwreal(x.^2, 200)
b = uwreal(x.^2, 200)
c = uwreal(x.^4, 200)
c = uwreal(x.^4, 200)
```
```
Correlations between variables are treated consistently. This requires that each variable that is defined with `uwreal` contains an ensemble `ID` tag. In the previous examples `a` has been measured on ensemble `ID` 1, while both `b` and `c` have been measured on ensemble `ID` 200. This will treat the measurements in `b` and `c` as statistically correlated, while `a` will be uncorrelated with both `b` and `c`.
Correlations between variables are treated consistently. Each call to `uwreal` contains an ensemble `ID` tag. In the previous examples `a` has been measured on ensemble `ID` 1, while both `b` and `c` have been measured on ensemble `ID` 200. Monte Carlo measurements with the same `ID` tag are assumed to be measured on the same configurations: `b` and `c` above are statistically correlated, while `a` will be uncorrelated with both `b` and `c`.
One can perform operations with `uwreal` variables as if there were normal floats
One can perform operations with `uwreal` variables as if there were normal floats
Error propagation is done automatically, and correlations are consistently taken into account. Note that once complex derived `uwreal` variables are defined, as for example `d` above, they will in general depend on more than one ensemble `ID`.
Error propagation is done automatically, and correlations are consistently taken into account. Note that once complex derived `uwreal` variables are defined, as for example `d` above, they will in general depend on more than one ensemble `ID`.
In order to perform the error analysis of a variable, one should use the `uwerr` function
The message "(Error not available... maybe run uwerr)" in all the above statements suggests that newly defined `uwreal` variables do not know their error yet. In order to perform the error analysis of a variable, one should use
the `uwerr` function
```@repl gs
```@repl gs
uwerr(d);
uwerr(d);
println("d: ", d)
println("d: ", d)
```
```
One can get detailed information on the error of a variables
A call to `uwerr` will apply the ``\Gamma``-method to MC ensembles. `ADerrors.jl` can give detailed information on the error of a variable (this requires that `uwerr` has been called)
```@repl gs
```@repl gs
uwerr(d);
println("Details on variable d: ")
println("Details on variable d: ")
details(d)
details(d)
```
```
where we can clearly see that there are two ensembles contributing to the uncertainty in `d`. We recognize this as `d` being a function of both `a` (measured on ensemble 1), and `b` and `c`, measured on ensemble 200.
Here we can see that there are two ensembles contributing to the uncertainty in `d`. We recognize this as `d` being a function of both `a` (measured on ensemble 1), and `b` and `c`, measured on ensemble 200. Most of the error in `d` comes in fact from ensemble 200.
Note that one does not need to use `uwerr` on a variable unless one is interested in the error on that variable. For example
Note that one does not need to use `uwerr` on a variable unless one is interested in the error on that variable. For example, continuing with the previous example
```@repl gs
```@repl gs
global x = 1.0
global x = 1.0
for i in 1:100
for i in 1:100
...
@@ -53,13 +57,13 @@ uwerr(x)
...
@@ -53,13 +57,13 @@ uwerr(x)
print("Root of d*cos(x) - x:" )
print("Root of d*cos(x) - x:" )
details(x)
details(x)
```
```
determines the root on ``d\cos(x) - x`` by using Newton's method, and propagates the error on `d` into an error on the root. The error on `x` is only determined once, after the 100 iterations are over.
determines the root on ``d\cos(x) - x`` by using Newton's method, and propagates the error on `d` into an error on the root. The error on `x` is only determined once, after the 100 iterations are over. This feature is useful becase calls to `uwerr` are typically numerically expensive in comparison with mathematical operations of `uwreal`'s.
## Advanced topics
## Advanced topics
### Errors in fit parameters
### Errors in fit parameters
`ADerrors.jl` does not provide an interface to perform fits, but once the minima of the ``\chi^2`` has been found, it can help propagating the error from the data to the fit parameters.
`ADerrors.jl` does not provide an interface to perform fits, but once the minima of the ``\chi^2`` has been found, it can return the fit parameters as `uwreal` variables.
Here we are going to [repeat the example](https://github.com/JuliaNLSolvers/LsqFit.jl) of the package `LsqFit.jl` using `ADerrors.jl` for error propagation.
Here we are going to [repeat the example](https://github.com/JuliaNLSolvers/LsqFit.jl) of the package `LsqFit.jl` using `ADerrors.jl` for error propagation.
`ADerrors.jl` can deal with observables that are not measured on every configuration, and still take correctly the correlations/autocorrelations into account. Here we show an extreme case where one observable is measured only in the even configurations, while the oher is measured on the odd coficurations.
`ADerrors.jl` can deal with observables that are not measured on every configuration, and still take correctly the correlations/autocorrelations into account. Here we show an extreme case where one observable is measured only in the even configurations, while the other is measured on the odd configurations.
In this case `uwerr` complains because with the chosen window the variance for ensemble with ID 1001 is negative. Looking at the normalized autocorrelation function we can easily spot the problem
the variables `x2` and `x4` are normal `uwreal` variables, and one can work with them as with any other variable. Note however that the automatic determination of the summation window can fail in these cases, because the autocorrelation function with missing measurements can be very complicated. For example
```@repl gaps
rat = x2/x4 # This is perfectly fine
uwerr(rat); # Use automatic window: might fail!
```
In this case `uwerr` complains because with the automatically chosen window the variance for ensemble with ID 1001 is negative. Let's have a look at the normalized autocorrelation function
```@repl gaps
```@repl gaps
iw = window(rat, 1001)
iw = window(rat, 1001)
r = rho(rat, 1001);
r = rho(rat, 1001);
...
@@ -141,19 +159,19 @@ dr = drho(rat, 1001);
...
@@ -141,19 +159,19 @@ dr = drho(rat, 1001);
plot(collect(1:100),
plot(collect(1:100),
r[1:100],
r[1:100],
yerr = dr[1:100],
yerr = dr[1:100],
seriestype = :scatter, title = "Chosen Window: " * string(iw))
The normalized autocorrelation function is oscillating, and the chosen window is 2! We better fix the window to 50 for this case
The normalized autocorrelation function is oscillating, and the automatically chosen window is 2, producing a negative variance! We better fix the window to 50 for this case
```@repl gaps
```@repl gaps
wpm = Dict{Int64, Vector{Float64}}()
wpm = Dict{Int64, Vector{Float64}}()
wpm[1001] = [50.0, -1.0, -1.0, -1.0]
wpm[1001] = [50.0, -1.0, -1.0, -1.0]
uwerr(rat, wpm)
uwerr(rat, wpm) # repeat error analysis with our choice (window=50)
println("Ratio: ", rat)
println("Ratio: ", rat)
```
```
After a visual examination of the autocorrelation function and a reasonable choice of the summation window, the error in `rat` is determined.
Note however that this is very observable dependent
Note however that it is very difficult to anticipate which observables will have a strange autocorrelation function. For example the product of `x2` and `x4` shows a reasonably well behaved autocorrelation function.
```@repl gaps
```@repl gaps
prod = x2*x4
prod = x2*x4
uwerr(prod)
uwerr(prod)
...
@@ -163,11 +181,11 @@ dr = drho(prod, 1001);
...
@@ -163,11 +181,11 @@ dr = drho(prod, 1001);
plot(collect(1:2*iw),
plot(collect(1:2*iw),
r[1:2*iw],
r[1:2*iw],
yerr = dr[1:2*iw],
yerr = dr[1:2*iw],
seriestype = :scatter, title = "Chosen Window: " * string(iw))
In general error analysis with data with arbitrary gaps is possible, and fully supported in `ADerrors.jl`, but it can be tricky, and certaiinly requires to examine the data carefully.
In general error analysis with data with arbitrary gaps is possible, and fully supported in `ADerrors.jl`. However in cases where many gaps are present, data analysis can be tricky, and certainly requires some care.
