Skip to content

Commit

Permalink
update tutorial
Browse files Browse the repository at this point in the history
  • Loading branch information
Cornelius-G committed Oct 27, 2023
1 parent a2dc5d7 commit 5ac79b5
Show file tree
Hide file tree
Showing 3 changed files with 33 additions and 240 deletions.
52 changes: 26 additions & 26 deletions docs/src/tutorial.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,12 @@ For our example, we consider two parameters with the names `C1` and `C2`.
For `C1` we choose a uniform (flat) prior in the range (-3, 3).
For `C2` we choose a gaussian prior with μ=0 and σ=0.5.

```julia
````julia
parameters = BAT.NamedTupleDist(
C1 = -3..3, # short for: Uniform(-3, 3)
C2 = Normal(0, 0.5) # Normal distribution
)
```
````

A parameter can be fixed (and therefore excluded from the fit) by setting its
prior to a certain value, e.g.: `C2 = 0`.
Expand All @@ -51,7 +51,7 @@ the parameters and that internally calls a more complex function and passes the
corresponding arguments. In this example, the function `xsec2` calls the
function `myfunc` and passes further arguments (`coeffs`).

```julia
````julia
function xsec1(params)
c = [20.12, 5.56, 325.556]
return c[1] * params.C1 + c[2] * params.C1 * params.C2+ c[3] * params.C2
Expand All @@ -65,13 +65,13 @@ end
function myfunc(params, c)
return c[1] * params.C1 + c[2] * params.C1 * params.C2+ c[3] * params.C2
end
```
````

If your observable is a distribution, you can define a vector of functions
with one function for each bin of the distribution.
(You could also treat each bin as a separate observable as shown above.)

```julia
````julia
function diff_xsec_bin1(params)
coeffs = [2.2, 5.5, 6.6]
return myfunc(params, coeffs)
Expand All @@ -88,7 +88,7 @@ function diff_xsec_bin3(params)
end

diff_xsec = [diff_xsec_bin1, diff_xsec_bin2, diff_xsec_bin3]
```
````

Note: Another way to define a vector of functions for the bins of distributions
is shown [here](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/advanced_tutorial/#Creating-a-vector-of-functions-for-distributions-1)
Expand Down Expand Up @@ -117,7 +117,7 @@ A `Measurement` can be excluded from the model by setting the switch `active=fal
For a `MeasurementDistribution`, the keyword `active` accepts `true` or `false`
to (de)activate the whole distribution or a vector of booleans for (de)activating only certain bins.

```julia
````julia
measurements = (
Meas1 = Measurement(xsec1, 21.6, uncertainties = (stat=0.8, syst=1.8, another_unc=2.3),
active=true), # `active = false`: exclude measurement from fit (default: active = true)
Expand All @@ -129,7 +129,7 @@ measurements = (
uncertainties = (stat = [0.7, 1.1, 1.2], syst= [0.7, 0.8, 1.3], another_unc = [1.0, 1.2, 1.9]),
active=[true, false, true]), # `active = false`: exclude all bins from fit, `active = [true, true, false]`: exclude only third bin from fit
)
```
````

Further information on the constructors see the API documentation of [`Measurement`](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/api/#EFTfitter.Measurement)
and [`MeasurementDistribution`](https://tudo-physik-e4.github.io/EFTfitter.jl/dev/api/#EFTfitter.MeasurementDistribution).
Expand Down Expand Up @@ -162,7 +162,7 @@ the Xth bin of the distribution can be accessed.
Note: This function is evaluated from top to bottom, so if you overwrite a
specific correlation value, the last value entered will be used.

```julia
````julia
dist_corr = [1.0 0.5 0.0;
0.5 1.0 0.0;
0.0 0.0 1.0]
Expand All @@ -173,9 +173,9 @@ another_corr_matrix = to_correlation_matrix(measurements,
(:MeasDist, :MeasDist, dist_corr), # correlate the bins of :MeasDist according to the matrix dist_corr
(:MeasDist_bin2, :MeasDist_bin3, 0.3), # correlate bin2 of :MeasDist with bin3 with 0.3 (overwrites the corresponding element set in the previous line, but ignored in fit since MeasDist_bin2 is inactive)
)
```
````

```julia
````julia
correlations = (
stat = NoCorrelation(active=true), # will use the identity matrix of the correct size

Expand All @@ -187,37 +187,37 @@ correlations = (

another_unc = Correlation(another_corr_matrix, active=true)
)
```
````

## File "runTutorial.jl"
Here, we create the `EFTfitterModel` from our inputs and run the actual analysis.

First, we need to setup EFTfitter, BAT and some other Julia packages:

```julia
````julia
using EFTfitter
using BAT # for sampling
using IntervalSets # for specifying the prior
using Distributions # for specifying the prior
using Plots # for plotting
```
````

We can then build the `EFTfitterModel` which combines all our inputs into
one object that is then used to perform the analysis on.

```julia
````julia
model = EFTfitterModel(parameters, measurements, correlations)
```
````

To sample the posterior distribution, we specify that our `EFTfitterModel`
should be used and then setup BAT.jl to sample the EFTfitter likelihood.

```julia
````julia
posterior = PosteriorMeasure(model)

algorithm = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)
samples = bat_sample(posterior, algorithm).result;
```
````

For further information on settings & algorithms when sampling with BAT.jl
see the BAT.jl [tutorial](https://bat.github.io/BAT.jl/dev/tutorial/#Parameter-Space-Exploration-via-MCMC)
Expand All @@ -226,7 +226,7 @@ and [documentation](https://bat.github.io/BAT.jl/dev/stable_api/#BAT.bat_sample)
We can then inspect the results of the sampling using BAT.jl's `SampledDensity`,
giving a summary of the sampling and the results of the model parameters.

```julia
````julia
sd = SampledDensity(posterior, samples)
display(sd)
```
Expand Down Expand Up @@ -258,16 +258,16 @@ cov ╲ │ C1 C2
───────┼─────────────────────────
C1 │ 0.122271 -0.0070394
C2 │ -0.0070394 0.000442548
```
````

Information about the smallest 1d intervals containing p% proability can be
obtained using the `get_smallest_interval_edges` function:

```julia
````julia
intervals_1d_C1 = get_smallest_interval_edges(samples, :C1, 0.9, bins=200, atol=0.1)
println("lower interval edges: $(intervals_1d_C1.lower)")
println("upper interval edges: $(intervals_1d_C1.upper)")
```
````

The keyword `atol` controls the absolute tolerance for which intervals are joined
together when they are seperated less than this value. This is particularly useful
Expand All @@ -276,21 +276,21 @@ when a large number of bins is used.
Of course, plotting the resulting posterior distributions is also simple
using Plots.jl and the BAT.jl plotting recipes:

```julia
````julia
p = plot(samples)
savefig(p, "plot.pdf")
```
````

![example plot](plots/plot.png)

For information on how to customize plots of the samples, please see the BAT.jl
[plotting documentation](https://bat.github.io/BAT.jl/dev/plotting/) and
[examples](https://github.com/bat/BAT.jl/blob/master/examples/dev-internal/plotting_examples.jl).

```julia
````julia
p = plot(samples, 0.9)
savefig(p, "plot_1d.pdf")
```
````

![example plot](plots/interval_plot_1.png)

Expand Down
Loading

0 comments on commit 5ac79b5

Please sign in to comment.