diff --git a/docs/source/design_example/animation.gif b/docs/source/design_example/animation.gif
new file mode 100644
index 0000000..1e5bac9
Binary files /dev/null and b/docs/source/design_example/animation.gif differ
diff --git a/docs/source/design_example/best_0.svg b/docs/source/design_example/best_0.svg
new file mode 100644
index 0000000..82fc293
--- /dev/null
+++ b/docs/source/design_example/best_0.svg
@@ -0,0 +1,889 @@
+
+
+
diff --git a/docs/source/design_example/best_1.svg b/docs/source/design_example/best_1.svg
new file mode 100644
index 0000000..bd3c252
--- /dev/null
+++ b/docs/source/design_example/best_1.svg
@@ -0,0 +1,884 @@
+
+
+
diff --git a/docs/source/design_example/loss.svg b/docs/source/design_example/loss.svg
new file mode 100644
index 0000000..6366c98
--- /dev/null
+++ b/docs/source/design_example/loss.svg
@@ -0,0 +1,873 @@
+
+
+
diff --git a/docs/source/design_example/parameters.svg b/docs/source/design_example/parameters.svg
new file mode 100644
index 0000000..4751e6f
--- /dev/null
+++ b/docs/source/design_example/parameters.svg
@@ -0,0 +1,1268 @@
+
+
+
diff --git a/docs/source/simple_stochastic_example/best_0.svg b/docs/source/simple_stochastic_example/best_0.svg
new file mode 100644
index 0000000..30be88c
--- /dev/null
+++ b/docs/source/simple_stochastic_example/best_0.svg
@@ -0,0 +1,1100 @@
+
+
+
diff --git a/docs/source/simple_stochastic_example/best_1.svg b/docs/source/simple_stochastic_example/best_1.svg
new file mode 100644
index 0000000..f40d91c
--- /dev/null
+++ b/docs/source/simple_stochastic_example/best_1.svg
@@ -0,0 +1,1606 @@
+
+
+
diff --git a/docs/source/simple_stochastic_example/best_2.svg b/docs/source/simple_stochastic_example/best_2.svg
new file mode 100644
index 0000000..abb58f1
--- /dev/null
+++ b/docs/source/simple_stochastic_example/best_2.svg
@@ -0,0 +1,1404 @@
+
+
+
diff --git a/examples/sample_curve_design/config.yaml b/examples/sample_curve_design/config.yaml
new file mode 100644
index 0000000..65aa7fb
--- /dev/null
+++ b/examples/sample_curve_design/config.yaml
@@ -0,0 +1,32 @@
+
+iters: 10
+
+optimiser: botorch
+
+
+parameters:
+ a: [0, -4, 4]
+ b: [0, -4, 4]
+ c: [0, -4, 4]
+
+
+
+objective:
+ name: design
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: abs( * exp(x) + * x**2 + * sin(x))
+ parametric: x
+ bounds: [-2, 2]
+ points: 100
+ targets:
+ 'minimum_value':
+ quantity: max
+ prediction: ['case_1']
+ negate: False
+ 'maximum_area':
+ quantity: integral
+ prediction: ['case_1']
+ negate: True
diff --git a/examples/sample_curve_design/description.md b/examples/sample_curve_design/description.md
new file mode 100644
index 0000000..c79cc08
--- /dev/null
+++ b/examples/sample_curve_design/description.md
@@ -0,0 +1,85 @@
+## sample_curve_design example
+
+A simple design objective example is included to demonstrate how to use `piglot` in design problems.
+
+In this problem, we aim to minimise the minimum value and maximimise the area below a nonlinear function given by the expression $f(x) = |a \exp(x) + bx^2+c\sin(x)|$.
+
+In this case, there is no reference response, and several values of the parameters `a`, `b` and `c` are evaluated for the target objectives.
+
+We run 10 iterations using the `botorch` optimiser (our interface for Bayesian optimisation), and set the parameters for optimisation (`a`, `b` and `c`) with bounds `[0,4]` and initial value 1.
+The notation ``, `` and `` indicates that these parameters should be optimised.
+We also define a parameterisation using the variable $x$, where we sample the function between `[-2,2]` with 100 points.
+
+### Using configuration files
+The configuration file (`examples/sample_curve_design/config.yaml`) for this example is:
+```yaml
+iters: 10
+
+optimiser: botorch
+
+parameters:
+ a: [0, -4, 4]
+ b: [0, -4, 4]
+ c: [0, -4, 4]
+
+objective:
+ name: design
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: abs( * exp(x) + * x**2 + * sin(x))
+ parametric: x
+ bounds: [-2, 2]
+ points: 100
+ targets:
+ 'minimum_value':
+ quantity: max
+ prediction: ['case_1']
+ negate: False
+ 'maximum_area':
+ quantity: integral
+ prediction: ['case_1']
+ negate: True
+```
+The generated response with the label `case_1` is optimised having two target objectives: (i) minimisation of the minimum function value `minimum_value`, and (ii) maximization of the area/integral below the function `maximum_area`.
+
+To run this example, open a terminal inside the `piglot` repository, enter the `examples/sample_curve_design` directory and run piglot with the given configuration file
+```bash
+cd examples/sample_curve_design
+piglot config.yaml
+```
+You should see an output similar to
+```
+BoTorch: 100%|███████████████████████████████████████| 10/10 [00:00<00:00, 10.92it/s, Loss: -4.3239e+00]
+Completed 10 iterations in 0.91609s
+Best loss: -4.32393814e+00
+Best parameters
+- a: -2.433969
+- b: -4.000000
+- c: 4.000000
+```
+In addition to these outputs, `piglot` creates an output directory, with the same name of the configuration file (minus the extension), where it stores the optimisation data.
+
+To visualise the optimisation results, use the `piglot-plot` utility.
+In the same directory, run (this may take some time)
+```bash
+piglot-plot animation config.yaml
+```
+This generates an animation for all the function evaluations that have been made throughout the optimisation procedure for the two target objectives.
+You can find the `.gif` file(s) inside the output directory, which should give something like (for the `minimum_value` target):
+![Best case plot](../../docs/source/design_example/animation.gif)
+
+Now, try running
+```bash
+piglot-plot parameters config.yaml
+```
+This will plot the evaluated parameters during the optimisation procedure:
+![Best case plot](../../docs/source/design_example/parameters.svg)
+
+To see the convergence history of the best loss function value, run
+```bash
+piglot-plot history config.yaml --best
+```
+which will generate:
+![Best case plot](../../docs/source/design_example/loss.svg)
\ No newline at end of file
diff --git a/examples/sample_curve_fitting/description.md b/examples/sample_curve_fitting/description.md
new file mode 100644
index 0000000..b6e01e4
--- /dev/null
+++ b/examples/sample_curve_fitting/description.md
@@ -0,0 +1,141 @@
+## sample_curve_fitting example
+
+A simple analytical curve fitting problem is included to demonstrate how to use `piglot`.
+In this case, we aim to fit a quadratic expression of the type $f(x) = a x^2$, using as a reference, a numerically generated reference from the expression $f(x) = 2 x^2$ (provided in the `examples/sample_curve_fitting/reference_curve.txt` file).
+We want to find the value for $a$ that better fits our reference (it should be 2).
+
+We run 10 iterations using the `botorch` optimiser (our interface for Bayesian optimisation), and set the parameter `a` for optimisation with bounds `[0,4]` and initial value 1.
+Our optimisation objective is the fitting of an analytical curve, with the expression ` * x ** 2`.
+The notation `` indicates that this parameter should be optimised.
+We also define a parameterisation using the variable $x$, where we sample the function between `[-5,5]` with 100 points.
+
+### Using configuration files
+The configuration file (`examples/sample_curve_fitting/config.yaml`) for this example is:
+```yaml
+iters: 10
+
+optimiser: botorch
+
+parameters:
+ a: [1, 0, 4]
+
+objective:
+ name: fitting
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ references:
+ 'reference_curve.txt':
+ prediction: ['case_1']
+```
+The generated response with the label `case_1` is compared with our reference response, given from the file `reference_curve.txt`
+
+To run this example, open a terminal inside the `piglot` repository, enter the `examples/sample_curve_fitting` directory and run piglot with the given configuration file
+```bash
+cd examples/sample_curve_fitting
+piglot config.yaml
+```
+You should see an output similar to
+```
+BoTorch: 100%|██████████████████████████████████████████████████████| 10/10 [00:00<00:00, 17.66it/s, Loss: 8.8505e-08]
+Completed 10 iterations in 0.56614s
+Best loss: 8.85050592e-08
+Best parameters
+- a: 1.999508
+```
+As you can see, piglot correctly identifies the `a` parameter close to the expected value of 2, and the error of the fitting is in the order of $10^{-8}$.
+In addition to these outputs, `piglot` creates an output directory, with the same name of the configuration file (minus the extension), where it stores the optimisation data.
+
+To visualise the optimisation results, use the `piglot-plot` utility.
+In the same directory, run
+```bash
+piglot-plot best config.yaml
+```
+Which will display the best observed value for the optimisation problem.
+You should see the following output in the terminal
+```
+Best run:
+Start Time /s 0.587397
+Run Time /s 0.004439
+a 1.999508
+Name: 18, dtype: object
+Hash: 2313718f75bc0445aa71df7d6d4e50ba82ad593d65f3762efdcbed01af338e30
+Objective: 8.85050592e-08
+```
+The script will also plot the best observed response, and its comparison with the reference response:
+![Best case plot](../../docs/source/simple_example/best.svg)
+
+Now, try running (this may take some time)
+```bash
+piglot-plot animation config.yaml
+```
+This generates an animation for all the function evaluations that have been made throughout the optimisation procedure.
+You can find the `.gif` file(s) inside the output directory, which should give something like:
+![Best case plot](../../docs/source/simple_example/animation.gif)
+
+
+### Using Python scripts
+
+Another way of using `piglot` is via its package and Python modules.
+This approach may offer increase flexibility in the setup of the optimisation problem, at the cost of increased complexity and verbosity.
+A sample script equivalent to the configuration file for the problem described in [the previous section](#using-configuration-files) is provided in `examples/sample_curve_fitting/config.py`, given by:
+```python
+import os
+import shutil
+from piglot.parameter import ParameterSet
+from piglot.solver.solver import Case
+from piglot.solver.curve.solver import CurveSolver
+from piglot.solver.curve.fields import CurveInputData, Curve
+from piglot.objectives.fitting import Reference, MSE
+from piglot.objectives.fitting import FittingObjective, FittingSolver
+from piglot.optimisers.botorch.bayes import BayesianBoTorch
+
+# Set up output and temporary directories
+output_dir = 'config'
+tmp_dir = os.path.join(output_dir, 'tmp')
+if os.path.isdir(output_dir):
+ shutil.rmtree(output_dir)
+os.makedirs(output_dir, exist_ok=True)
+
+# Set up optimisation parameters
+parameters = ParameterSet()
+parameters.add('a', 1.0, 0.0, 4.0)
+
+# Set up the reference
+reference = Reference('reference_curve.txt', ['case_1'], output_dir)
+
+# Set up the solver to use
+input_data = CurveInputData('case_1', ' * x ** 2', 'x', (-5.0, 5.0), 100)
+case_1 = Case(input_data, {'case_1': Curve()})
+solver = CurveSolver([case_1], parameters, output_dir, tmp_dir=tmp_dir)
+
+# Set up the fitting objective
+references = {reference: ['case_1']}
+fitting_solver = FittingSolver(solver, references)
+objective = FittingObjective(parameters, fitting_solver, output_dir, MSE())
+
+# Set up the optimiser and run optimisation
+optimiser = BayesianBoTorch(objective)
+value, params = optimiser.optimise(10, parameters, output_dir)
+print(f"Optimal value: {value}")
+print(f"Optimal parameters: {params}")
+```
+Run with
+```bash
+python config.py
+```
+Example output
+```
+BoTorch: 100%|██████████████████████████████████████████████████████| 10/10 [00:00<00:00, 16.75it/s, Loss: 8.9167e-08]
+Completed 10 iterations in 0.59692s
+Best loss: 8.91673999e-08
+Best parameters
+- a: 1.999506
+Optimal value: 8.916739991036405e-08
+Optimal parameters: [1.99950592]
+```
diff --git a/examples/sample_curve_fitting_composite/config.yaml b/examples/sample_curve_fitting_composite/config.yaml
new file mode 100644
index 0000000..031a9b5
--- /dev/null
+++ b/examples/sample_curve_fitting_composite/config.yaml
@@ -0,0 +1,21 @@
+iters: 10
+
+optimiser: botorch
+
+parameters:
+ a: [1, 0, 4]
+
+objective:
+ name: fitting
+ composite: True
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ references:
+ 'reference_curve.txt':
+ prediction: ['case_1']
diff --git a/examples/sample_curve_fitting_composite/description.md b/examples/sample_curve_fitting_composite/description.md
new file mode 100644
index 0000000..eb4d948
--- /dev/null
+++ b/examples/sample_curve_fitting_composite/description.md
@@ -0,0 +1,89 @@
+## sample_curve_fitting_composite example
+
+A simple analytical curve fitting problem using a composite strategy is included to demonstrate how to use `piglot` in the composite setting.
+In this case, we aim to fit a quadratic expression of the type $f(x) = a x^2$, using as a reference, a numerically generated reference from the expression $f(x) = 2 x^2$ (provided in the `examples/sample_curve_fitting/reference_curve.txt` file).
+We want to find the value for $a$ that better fits our reference (it should be 2).
+
+We run 10 iterations using the `botorch` optimiser (our interface for Bayesian optimisation), and set the parameter `a` for optimisation with bounds `[0,4]` and initial value 1.
+Our optimisation objective is the fitting of an analytical curve, with the expression ` * x ** 2`.
+The notation `` indicates that this parameter should be optimised.
+We also define a parameterisation using the variable $x$, where we sample the function between `[-5,5]` with 100 points.
+
+The particularity of this example is that a composite strategy is used to fit the response.
+The advantages of this composite Bayesian optimisation are demonstrated in [Coelho et al.]([docs/source/simple_example/best.svg](https://dx.doi.org/10.2139/ssrn.4674421)) and are two-fold: (i) more accurate posteriors for the loss function and (ii) reduced loss of information during the computation of the reduction function.
+
+In short, in the composite setting the loss function $\mathcal{L}\left(\bm{\theta}\right)$ is written as
+$
+ \mathcal{L}\left(\bm{\theta}\right)
+ =
+ \hat{\mathcal{L}}\left(\bm{e}\left(\bm{\theta}\right)\right)
+ =
+ \dfrac{1}{N}
+ \sum_{i=1}^{N}
+ \left[e_i\left(\bm{\theta}\right)\right]^2,
+$
+where $\bm{e}$ is a vector containing the pointwise errors at every reference point, $\hat{\mathcal{L}}\left(\bullet\right)$ is the scalar reduction function applied (NMSE in this case) and $N$ is the number of reference points.
+Thus, the problem can be stated as the minimisation of a composite function $\hat{\mathcal{L}}\left(\bm{e}\left(\bm{\theta}\right)\right)$, where $\hat{\mathcal{L}}\left(\bullet\right)$ is known and $\bm{e}\left(\bm{\theta}\right)$ is unknown.
+
+Consider the minimisation of $\mathcal{L}\left(\bm{\theta}\right) = \hat{\mathcal{L}}\left(\bm{e}\left(\bm{\theta}\right)\right)$.
+Within this setting, we start by replacing the single-output Gaussian Process on the loss $\mathcal{L}\left(\bm{\theta}\right)$ with a multi-output Gaussian Process for the error of each reference point $e_i$, that is,
+$
+ e_i\left(\bm{\theta}\right)
+ \sim
+ \mathcal{GP}
+ \left(
+ \mu_i\left(\bm{\theta}\right),
+ k_i\left(\bm{\theta},\bm{\theta}'\right)
+ \right).
+$
+Naturally, this requires feeding the optimiser with the entire response instead of a single scalar value.
+At this stage, each GP is assumed independent and the correlations between the outputs are not considered; that is, each value $e_i$ is assumed as independent and uniquely defined by the set of parameters $\bm{\theta}$.
+
+The configuration file (`examples/sample_curve_fitting_composite/config.yaml`) for this example is:
+```yaml
+iters: 10
+
+optimiser: botorch
+
+parameters:
+ a: [1, 0, 4]
+
+objective:
+ name: fitting
+ composite: True
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ references:
+ 'reference_curve.txt':
+ prediction: ['case_1']
+```
+The composite strategy is activated by setting ```composite: True```.
+
+To run this example, open a terminal inside the `piglot` repository, enter the `examples/sample_curve_fitting_composite` directory and run piglot with the given configuration file
+```bash
+cd examples/sample_curve_fitting_composite
+piglot config.yaml
+```
+You should see an output similar to
+```
+BoTorch: 100%|████████████████████████████████████████| 10/10 [00:01<00:00, 7.94it/s, Loss: 5.6009e-08]
+Completed 10 iterations in 1s
+Best loss: 5.60089334e-08
+Best parameters
+- a: 1.999685
+```
+It is observed that piglot correctly identifies the `a` parameter close to the expected value of 2, and the error of the fitting is in the order of $10^{-8}$.
+As this example is quite simple, there is no great advantage of using the composite strategy, as the simple Bayesian optimisation is already able of finding accurate solutions within few function evaluations.
+
+
+To visualise the optimisation results, use the `piglot-plot` utility.
+In the same directory, run for instance
+```bash
+piglot-plot best config.yaml
+```
\ No newline at end of file
diff --git a/examples/sample_curve_fitting_composite/reference_curve.txt b/examples/sample_curve_fitting_composite/reference_curve.txt
new file mode 100644
index 0000000..3d109ab
--- /dev/null
+++ b/examples/sample_curve_fitting_composite/reference_curve.txt
@@ -0,0 +1,6 @@
+-4.0 32.0
+-2.4 11.52
+-0.7999999999999998 1.2799999999999994
+0.8000000000000007 1.2800000000000022
+2.4000000000000004 11.520000000000003
+4.0 32.0
diff --git a/examples/sample_curve_fitting_stochastic/config.yaml b/examples/sample_curve_fitting_stochastic/config.yaml
new file mode 100644
index 0000000..bbd9c6c
--- /dev/null
+++ b/examples/sample_curve_fitting_stochastic/config.yaml
@@ -0,0 +1,27 @@
+iters: 10
+
+optimiser: botorch
+
+parameters:
+ a: [1, 0, 4]
+
+objective:
+ name: fitting
+ stochastic: True
+ composite: False
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ 'case_2':
+ expression: 2* * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ references:
+ 'reference_curve.txt':
+ prediction: ['case_1', 'case_2']
\ No newline at end of file
diff --git a/examples/sample_curve_fitting_stochastic/description.md b/examples/sample_curve_fitting_stochastic/description.md
new file mode 100644
index 0000000..0614b42
--- /dev/null
+++ b/examples/sample_curve_fitting_stochastic/description.md
@@ -0,0 +1,87 @@
+## sample_curve_fitting_stochastic example
+
+A simple analytical curve fitting problem with noise in the input data is included to demonstrate how to use `piglot` with variance.
+In this case, we aim to fit two quadratic expressions of the type $f(x) = a x^2$ and $f(x) = 2a x^2$, using as a reference response, a numerically generated reference from the expression $f(x) = 2 x^2$ (provided in the `examples/sample_curve_fitting_stochastic/reference_curve.txt` file).
+We want to find the value for $a$ that better fits our reference. In this case, the optimal solution is no longer $a=2$, as two distinct functions are used to compute the generated response.
+
+We run 10 iterations using the `botorch` optimiser (our interface for Bayesian optimisation), and set the parameter `a` for optimisation with bounds `[0,4]` and initial value 1.
+The notation `` indicates that this parameter should be optimised.
+We also define a parameterisation using the variable $x$, where we sample the function between `[-5,5]` with 100 points.
+
+
+DESCRIBE STOCHASTIC
+
+The configuration file (`examples/sample_curve_fitting_stochastic/config.yaml`) for this example is:
+```yaml
+iters: 10
+
+optimiser: botorch
+
+parameters:
+ a: [1, 0, 4]
+
+objective:
+ name: fitting
+ stochastic: True
+ composite: False
+ solver:
+ name: curve
+ cases:
+ 'case_1':
+ expression: * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ 'case_2':
+ expression: 2* * x ** 2
+ parametric: x
+ bounds: [-5, 5]
+ points: 100
+ references:
+ 'reference_curve.txt':
+ prediction: ['case_1', 'case_2']
+```
+The stochastic strategy is activated by setting ```stochastic: True```, and by adding a new generated response with the label `case_2`, given by the expression $f(x) = 2a x^2$.
+
+To run this example, open a terminal inside the `piglot` repository, enter the `examples/sample_curve_fitting_stochastic` directory and run piglot with the given configuration file
+```bash
+cd examples/sample_curve_fitting_stochastic
+piglot config.yaml
+```
+You should see an output similar to
+```
+BoTorch: 100%|████████████████████████████████████████| 10/10 [00:00<00:00, 14.42it/s, Loss: 1.7315e-01]
+Completed 10 iterations in 0.69363s
+Best loss: 1.73146009e-01
+Best parameters
+- a: 1.198191
+```
+Piglot identifies the `a` parameter as 1.2, and the error of the fitting is in the order of $10^{-1}$.
+For this case, the use of the composite Bayesian strategy (as decribed in [here](examples/sample_curve_fitting_composite/description.md)) significantly improves the quality of the fitting.
+
+To visualise the optimisation results, use the `piglot-plot` utility.
+In the same directory, run
+```bash
+piglot-plot best config.yaml
+```
+Which will display the best observed value for the optimisation problem.
+You should see the following output in the terminal
+```
+Best run:
+Start Time /s 0.683277
+Run Time /s 0.014079
+Variance 0.005519
+a 1.198191
+Name: 18, dtype: object
+Hash: f07c094fdbbaa637387a31cdeeb946783f4b3aeefe99639995d7e6539cf48475
+Objective: 1.73146009e-01
+```
+The script plots the best observed responses, and its comparison with the reference response. Moreover, the average, the median, the standard deviation and the 95\% mean condidence intervals are also provided.
+![Best case plot](../../docs/source/simple_stochastic_example/best_0.svg)
+![Best case plot](../../docs/source/simple_stochastic_example/best_1.svg)
+![Best case plot](../../docs/source/simple_stochastic_example/best_2.svg)
+
+
+Try running the same example with the composite strategy (by simply setting ```composite: True```). With the composite Bayesian optimisation the error of the fitting is reduced to $10^{-5}$.
+
+
diff --git a/examples/sample_curve_fitting_stochastic/reference_curve.txt b/examples/sample_curve_fitting_stochastic/reference_curve.txt
new file mode 100644
index 0000000..3d109ab
--- /dev/null
+++ b/examples/sample_curve_fitting_stochastic/reference_curve.txt
@@ -0,0 +1,6 @@
+-4.0 32.0
+-2.4 11.52
+-0.7999999999999998 1.2799999999999994
+0.8000000000000007 1.2800000000000022
+2.4000000000000004 11.520000000000003
+4.0 32.0