diff --git a/docs/Project.toml b/docs/Project.toml
index 92da694..7be67d5 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,8 @@
 [deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SafetySignalDetection = "ad9a3248-762c-4ade-a843-8f185e2c18ff"
diff --git a/docs/src/introduction.md b/docs/src/introduction.md
new file mode 100644
index 0000000..3dec2a4
--- /dev/null
+++ b/docs/src/introduction.md
@@ -0,0 +1,110 @@
+# Introduction
+
+We will introduce the use of the `SafetySignalDetection.jl` package with a small example.
+
+## Data example
+
+We use a small example data set. 
+
+Note that:
+
+- The outcome `y` needs to be a bool vector. A different type of vector (e.g. integer) will fail later for the model fitting.
+- The column `time` needs to be a `Float64` (time until adverse event or until last treatment or follow up).
+- The column `trialindex` needs to be `Int64` (index of historical trials, starting from 1 and consecutively numbered).
+
+```julia
+using CSV
+using DataFrames
+using SafetySignalDetection
+
+module_dir = dirname(dirname(pathof(SafetySignalDetection)))
+csv_path = joinpath(module_dir, "test", "small_historic_dataset.csv")
+df = DataFrame(CSV.File(csv_path))
+df.y = map(x -> x == 1, df.y)
+if ("trial" ∈ names(df))
+    rename!(df, [:trial] .=> [:trialindex])
+end
+```
+
+## Meta-analytic prior sampling
+
+For both parameters $a$ and $b$, we can use beta distributions, specified via `Distributions.Beta()`.
+
+For the mean $a$, we can set the prior mean according to the background rate of the adverse event. Note that this is per unit of time. So depending how you receive the incidence rate in which units, you need to convert this back to the time unit you are using for the analysis.
+The overall sum of the two beta distribution's parameters can be kept constant at a low number.
+
+For the discount factor $b$, we can set the prior for it closer to 1 for more homogenous trials,
+and closer to 0 for more heterogenous trials.
+For both priors, the prior information should be low, i.e. the sum of the two parameters
+of the beta distribution should be not high, e.g. around 10.
+
+Note:
+
+- The number of samples (e.g. 10,000) and the number of threads (1) need to be passed as integers. Otherwise the call will fail.
+- We don't need to discard a burn in period now because the NUTS sampler is used internally and is discarding a burn in automatically already.
+
+```julia
+using Distributions
+
+prior_a = Beta(1 / 3, 1 / 3)
+prior_b = Beta(5, 5)
+
+map_samples = meta_analytic_samples(df, prior_a, prior_b, 10_000, 1)
+```
+
+The resulting $\pi^{*}$ samples in `map_samples` are from the meta-analytic prior for the adverse event rate per unit of time in the control arm in the ongoing blinded trial.
+
+## Closed form approximation
+
+Now we give the $\pi^{*}$ samples to the `fit_beta_mixture()` function, together with the 
+number of components. Usually 2 components will be sufficient to achieve a good approximation.
+Internally, the classic Expectation Maximization (EM) algorithm is used.
+
+We can check the approximation graphically by overlaying the resulting probability density function with the
+samples histogram.
+
+```julia
+using Plots
+
+map_approx = fit_beta_mixture(map_samples, 2)
+
+stephist(map_samples, label = "prior samples", norm = :pdf)
+xvals = range(minimum(map_samples), maximum(map_samples); length = 100)
+plot!(xvals, pdf(map_approx, xvals), label = "approximate prior")
+```
+
+## Blinded trial analysis
+
+Now we can proceed to analyzing the blinded clinical trial.
+
+The important new ingredients here are:
+
+- The prior on the experimental adverse event rate. This can be uninformative, or weakly informative using a diluted background rate distribution.
+- The proportion of patients in the experimental treatment arm.
+
+```julia
+prior_exp = Beta(1, 1)
+prior_ctrl = map_approx
+exp_proportion = 0.5
+
+csv_current_path = joinpath(module_dir, "docs", "src", "small_current_dataset.csv")
+df_current = DataFrame(CSV.File(csv_current_path))
+df_current.y = map(x -> x == 1, df_current.y)
+
+result = blinded_analysis_samples(df_current, prior_exp, prior_ctrl, exp_proportion, 10_000, 1)
+```
+
+## Summary
+
+We can now e.g. look at the posterior probability that the rate in the experimental arm, $\pi_E$, is larger than the event rate in the control arm, $\pi_C$.
+
+```julia
+mean(result[!, :pi_exp] .> result[!, :pi_ctrl])
+```
+
+We can also create corresponding plots.
+
+```julia
+stephist(result[!, :pi_exp], label = "experimental arm", norm = :pdf)
+stephist!(result[!, :pi_ctrl], label = "control arm", norm = :pdf)
+```
diff --git a/docs/src/small_current_dataset.csv b/docs/src/small_current_dataset.csv
new file mode 100644
index 0000000..8acba6f
--- /dev/null
+++ b/docs/src/small_current_dataset.csv
@@ -0,0 +1,101 @@
+"","y","time","tmt"
+"1",0,0.508060037142824,2
+"2",0,0.5208363096442,2
+"3",0,0.556331578741567,2
+"4",0,0.563791380921675,1
+"5",0,0.571754800768225,1
+"6",0,0.590052150745623,1
+"7",0,0.591040918682293,2
+"8",0,0.595221235355885,2
+"9",0,0.597784269735739,1
+"10",0,0.607071721277385,1
+"11",0,0.611330809521675,2
+"12",1,0.613947846980951,2
+"13",0,0.618876835716143,2
+"14",1,0.640354126490284,2
+"15",0,0.645504716144507,2
+"16",0,0.648776558221881,2
+"17",0,0.653263780100834,2
+"18",0,0.664869826653364,2
+"19",0,0.671171271026457,2
+"20",0,0.674994900783847,1
+"21",1,0.679323147659513,2
+"22",0,0.681681905904048,2
+"23",0,0.693983536523541,1
+"24",0,0.713401719711471,1
+"25",0,0.715330274447549,1
+"26",1,0.715554267597332,2
+"27",0,0.724320719842258,1
+"28",0,0.743751407473829,2
+"29",0,0.744621157756557,2
+"30",0,0.749986377628643,2
+"31",0,0.756438894458238,2
+"32",1,0.769038580022455,1
+"33",0,0.781380031033633,2
+"34",0,0.782826549365607,2
+"35",0,0.78411009277818,2
+"36",1,0.791492904824573,2
+"37",1,0.801637802299991,1
+"38",0,0.811791839104605,2
+"39",1,0.837362781844232,1
+"40",0,0.841242652393076,1
+"41",0,0.862928944447668,2
+"42",0,0.869636503880502,1
+"43",0,0.886621511367525,2
+"44",0,0.889937450386028,2
+"45",0,0.89081983407345,2
+"46",0,0.894065672679249,2
+"47",0,0.895545855571809,1
+"48",0,0.909240301239177,2
+"49",0,0.919524707524235,2
+"50",0,0.924650129098574,2
+"51",0,0.945304040451132,2
+"52",0,0.946978089384413,1
+"53",0,0.951859002214535,1
+"54",0,0.957347076003118,2
+"55",0,0.960843806894175,1
+"56",1,0.965779375707409,2
+"57",0,0.966733577783565,1
+"58",0,0.96968774609861,2
+"59",0,0.972657807600261,2
+"60",1,0.979699684335007,1
+"61",0,0.979834431354125,2
+"62",0,0.990429732491942,2
+"63",0,0.993043940746905,2
+"64",0,0.9954099067624,2
+"65",0,0.997159593415083,2
+"66",0,1.00392896736284,1
+"67",0,1.01355110456085,2
+"68",0,1.04904199831044,2
+"69",0,1.06484850264858,1
+"70",0,1.06794773782012,2
+"71",1,1.0712293948771,2
+"72",0,1.08236479130322,1
+"73",0,1.0874704878659,2
+"74",0,1.08815656886374,1
+"75",0,1.09446561418411,2
+"76",1,1.11639969562923,2
+"77",1,1.13149712565995,1
+"78",1,1.15531456740156,2
+"79",0,1.17462593214085,2
+"80",1,1.18276276532017,2
+"81",0,1.18776627686728,1
+"82",0,1.20208448878763,2
+"83",1,1.21019494540132,2
+"84",0,1.21124169832635,1
+"85",0,1.2113759552136,1
+"86",0,1.21897941682307,1
+"87",0,1.22070786681776,2
+"88",0,1.24966304313974,1
+"89",0,1.2527988999829,1
+"90",0,1.25518082340978,2
+"91",0,1.25891242750253,2
+"92",0,1.26357294799558,2
+"93",0,1.26469715973161,1
+"94",0,1.27593270482697,2
+"95",0,1.28273087209758,2
+"96",0,1.29175769180613,1
+"97",0,1.31240415124692,2
+"98",1,1.31259288685395,1
+"99",0,1.33657055292152,2
+"100",1,1.34760702875274,2
diff --git a/test/large_historic_dataset.csv b/test/large_historic_dataset.csv
index 9f65459..428483f 100644
--- a/test/large_historic_dataset.csv
+++ b/test/large_historic_dataset.csv
@@ -1,4 +1,4 @@
-"","y","trial","time"
+"","y","trialindex","time"
 "1",0,1,0.508060037142824
 "2",0,7,0.5208363096442
 "3",0,2,0.556331578741567
diff --git a/test/small_historic_dataset.csv b/test/small_historic_dataset.csv
index bf51bee..7517dbc 100644
--- a/test/small_historic_dataset.csv
+++ b/test/small_historic_dataset.csv
@@ -1,4 +1,4 @@
-"","y","trial","time"
+"","y","trialindex","time"
 "1",0,5,0.508060037142824
 "2",0,4,0.5208363096442
 "3",0,2,0.556331578741567
diff --git a/test/test_meta_analysis.jl b/test/test_meta_analysis.jl
index 11d55a7..56d0c53 100644
--- a/test/test_meta_analysis.jl
+++ b/test/test_meta_analysis.jl
@@ -32,7 +32,7 @@ end
     rng = StableRNG(123)
     map_small = sample(
         rng,
-        meta_analysis_model(df_small.y, df_small.time, df_small.trial, 
+        meta_analysis_model(df_small.y, df_small.time, df_small.trialindex, 
                             prior_a, prior_b), 
         NUTS(0.65),
         10_000
@@ -54,7 +54,7 @@ end
     rng = StableRNG(123)
     map_large = sample(
         rng,
-        meta_analysis_model(df_large.y, df_large.time, df_large.trial, 
+        meta_analysis_model(df_large.y, df_large.time, df_large.trialindex, 
                             prior_a, prior_b), 
         NUTS(0.65), 
         10_000