Exercise-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BAT.jl Tutorial - Poisson Counting Experiment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: This tutorial is written in the style of a cloze text. Parts that have to be filled out by you are marked `#=FILLME=#`.\n",
    "\n",
    "The [the tutorial section](https://bat.github.io/BAT.jl/stable/tutorial/) of the BAT.jl documentation should prove helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BAT\n",
    "using Distributions \n",
    "using IntervalSets\n",
    "using ValueShapes\n",
    "using Plots\n",
    "using ArraysOfArrays\n",
    "using StatsBase \n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Situation\n",
    "We want to determine the properties of a radioactive singal source in the presence of background from natural sources of radioactivity.\n",
    "We assume to have one signal source $S$ and only one source of background $B$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Background only measurement\n",
    "We start by using our detector without the signal source installed.\n",
    "This measurement yields a number of $k_B=10$ counts.  \n",
    "Using this measurement we want to gain information about the event rate of the natural radioactive background.\n",
    "\n",
    "### Task: \n",
    "Perform a Bayesian analysis to estimate the event rate of the natural background $\\lambda_b$ using a Poisson model. \n",
    "\n",
    "Start by defining the log-likelihood function, using the function `logpdf()` and type `Poisson` provided by the package [Distributions.jl](https://juliastats.github.io/Distributions.jl/latest/univariate/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of observed background events\n",
    "kb = 10\n",
    "\n",
    "likelihood_B = let k = kb\n",
    "    params -> begin\n",
    "        return logpdf(Poisson(params.λb), k) # poisson log-likelihood\n",
    "    end\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the Prior with help of the `NamedTupleDist()` function. Use a flat prior between 0 and 30:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_B = NamedTupleDist(\n",
    "    λb = #=FILLME=#\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then use the likelihood and the prior to define the `PosteriorDensity()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_B = PosteriorDensity(#=FILLME=#); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the settings for the sampling. Choose `MetropolisHastings()` as the algorithm and set the number of chains and samples (4 chains with 10^5 samples each should be fine):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "algorithm = #=FILLME=#)\n",
    "nchains = #=FILLME=#)\n",
    "nsamples = =#=FILLME=#)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the function `bat_sample` to sample the posterior using `nchains` MCMC chains, generating `nsamples` per chain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_B, stats_B = bat_sample(posterior_B, (nsamples, nchains), algorithm)\n",
    "stats_B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the resulting disribution for the background rate using `plot()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot(posterior_B, samples_B, :λb)\n",
    "plot!(prior_B, :λb, xlabel = \"\\$\\\\lambda_b\\$\", ylabel = \"\\$P(\\\\lambda_b)\\$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print some statistics of the samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"Mode: $(stats_B.mode)\")\n",
    "println(\"Mean: $(stats_B.mean)\")\n",
    "println(\"Covariance: $(stats_B.cov)\")\n",
    "println(\"Standard Deviation: $(diag(sqrt(stats_B.cov[:, :])))\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Second background only measurement\n",
    "A second measurement of the natural background yields a number of $k_B=8$ counts.  \n",
    "Therefore, we want to update our estimation of the background rate using this new knowledge together with the pervious results.\n",
    "### Task:\n",
    "Perform an analysis of the new measurement similar to the first one.   \n",
    "Use the posterior distribution of the previous background measurement as a prior for this analysis.  \n",
    "This can be done by using a [StatsBase histogram](http://juliastats.github.io/StatsBase.jl/latest/empirical/#Histograms-1) of the samples: `fit(Histogram, flatview(samples)[i,:], weights, nbins)`.  \n",
    "The weights of the samples need to be used with the `FrequencyWeights()` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a `StatsBase` histogram containing the previous posterior distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_hist_B1 = fit(Histogram, flatview(samples_B.params)[1, :], FrequencyWeights(samples_B.weight), nbins = 400);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histogram can be used as a prior by converting it into a univariate distribution using `BAT.HistogramAsUvDistribution(histogram)`.  \n",
    "Apart from that, you can proceed similarly to the first task."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the log-likelihood function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kb2 = 8\n",
    "\n",
    "likelihood_B2 = let k = kb2\n",
    "    params -> #=FILLME=#\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the prior (use the posterior of the previous task) and the posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_B2 = NamedTupleDist( #=FILLME=# )\n",
    "\n",
    "posterior_B2 = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_B2, stats_B2 = #=FILLME=#\n",
    "stats_B2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the  `plot(samples)` and `plot!(prior)` functions to visualize the posterior of the first analysis and the updated posterior together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(posterior_B2, #=FILLME=#)\n",
    "plot!(prior_B2, #=FILLME=#)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print some statistics of the samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"Mode: #=FILLME=#\")\n",
    "println(\"Mean: #=FILLME=#\")\n",
    "println(\"Covariance: #=FILLME=#\")\n",
    "println(\"Standard Deviation: #=FILLME=#\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Signal + Background\n",
    "Including the radioactive source to the experimental setup, we repeat the measurement and obtain $k_{S+B}=12$ counts.\n",
    "With this measurement and our prior knowledge about the background we are able to estimate the rate of the signal $\\lambda_s$.\n",
    "### Task\n",
    "Perform a third analysis using a poisson model with the combined singal + background rate.\n",
    "Use the known information about the background from the previous tasks as a prior and choose a suitable prior for the signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the likelihood for the signal + background model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of observed events\n",
    "kSB = 12\n",
    "\n",
    "likelihood_SB = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the prior for both the signal and backgound parameters.  \n",
    "Remember to use the knowledge from the previous tasks for the background. (hint: `BAT.HistogramAsUvDistribution()`)  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_B2 = #=FILLME=#\n",
    "B2 = #=FILLME=#\n",
    "\n",
    "prior_SB = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also define the posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior_SB = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate samples for the signal + background model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_SB = #=FILLME=#\n",
    "stats_SB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize an overview of the results for both prameters use `plot(samples)`.   \n",
    "(Hint: use the keyword `param_labels=[\"\\\\lambda_b\",\"\\\\lambda_s\"]` for correct axis labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(samples_SB, #=FILLME=#)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print some statistics of the samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"Mode: #=FILLME=#\")\n",
    "println(\"Mean: #=FILLME=#\")\n",
    "println(\"Covariance: #=FILLME=#\")\n",
    "println(\"Standard Deviation: #=FILLME=#\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Error propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we want to caluclate the cross section of the signal process.\n",
    "The rate of measured events in the detector of a couting experiment can be written as \n",
    "\n",
    "$\\frac{\\mathrm d N}{\\mathrm d t} = \\epsilon \\cdot σ \\cdot L$ ,\n",
    "\n",
    "with the Luminosity $L$ and the efficiency of the detector $\\epsilon$.   \n",
    "The signal cross section is therefore given as:\n",
    "### $σ_S = \\frac{λ_s}{\\epsilon \\cdot L}$  .\n",
    "For this experiment we assume a luminosity $L = 1.1$  (neglecting units).\n",
    "\n",
    "As a final result we want to obtain either a measurement or an upper limit on the signal cross section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4 a) Known efficiency with gaussian uncertainty\n",
    "The detector efficiency has been measured to be $\\epsilon = 0.1 \\pm 0.02$, assuming the uncertainties to follow a normal distribution.\n",
    "\n",
    "Calculate the signal cross section $σ_S$ using the equation above.\n",
    "\n",
    "Use the [Distrtibutions.jl](https://juliastats.github.io/Distributions.jl/latest/univariate/) package and `rand()` to obain a sample for $\\epsilon$.  \n",
    "In order to obtain unweigthed samples of $\\lambda_S$ by resampling the posterior samples, the function `bat_sample(samples, nsamples).samples` can be used.\n",
    "\n",
    "The function `broadcast()` or the `.`operator (e.g. `a .+ b`) might be useful for element-wise operation when handeling the samples.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the luminosity and the efficiency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = #=FILLME=#\n",
    "ϵ = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the efficiency.  \n",
    "Use a [StatsBase histogram](http://juliastats.github.io/StatsBase.jl/latest/empirical/#Histograms-1) to visualize the distribution.  \n",
    "(Hint: The plot recipes can also be used for the `StatsBase` histograms) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_ϵ = fit(Histogram, ϵ, nbins=200)\n",
    "#!EX hist_ϵ = fit(#=FILLME=#)\n",
    "\n",
    "plot(hist_ϵ, 1, seriestype=:smallest_intervals, xlabel=\"\\$\\\\epsilon\\$\", ylabel=\"\\$p(\\\\epsilon)\\$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get unweighted samples for the signal rate and calculate the cross section distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resampled_SB = bat_sample(samples_SB, nsamples).samples;\n",
    "λ_SB = valshape(prior_SB).(resampled_SB.params).λs\n",
    "σS = (λ_SB)./(ϵ*L);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the distribution of the signal cross section. (Hint: use again a histogram)  \n",
    "From the plot, determine the 95% upper limit on the cross section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_σ = #=FILLME=#\n",
    "plot(hist_σ, #=FILLME=#)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4 b) Binomial analysis of calibration measurement to determine efficiency \n",
    "We now want to perform the same analysis as in task 4 a) but for the case that the detector efficiency $\\epsilon$ is not yet known.   \n",
    "\n",
    " Instead, $\\epsilon$ is to be determined using a calibration measurement with a source for which the signal rate is known. Then it is possible to calculate the efficiency from the expected number of counts and the number of counts actually measured with the detector.   \n",
    "\n",
    "\n",
    "For our example, the number of expected events is assumed to be $N_\\text{expected} = 200$.  \n",
    "The detector measures only $N_\\text{measured} = 21$ events.  \n",
    "\n",
    "Task: Implement a binomial model using the `Binomial(n,p)` function of the [Distrtibutions.jl](https://juliastats.github.io/Distributions.jl/latest/univariate/) package and determine the distribution of the detector efficiency.  \n",
    "Afterwards, repeat the steps from 4 a) using the obtained distribution of the efficiency to calculate the signal cross section ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the binomial likelihood:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_expected = 200\n",
    "n_measured = 21\n",
    "\n",
    "likelihood_binomial = let #=FILLME=#\n",
    "    params -> #=FILLME=#\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the prior (flat) for the efficiency and define the posterior:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_binomial = #=FILLME=#\n",
    "\n",
    "posterior_binomial = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_binomial, stats_binomial = #=FILLME=#\n",
    "stats_binomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the distribution of the efficiency:  \n",
    "(What do you observe when comparing to the efficiency used in 4 a) ?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(samples_binomial, #=FILLME=#)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print some statistics of the samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"Mode: $(stats_binomial.mode)\")\n",
    "println(\"Mean: $(stats_binomial.mean)\")\n",
    "println(\"Covariance: $(stats_binomial.cov)\")\n",
    "println(\"Standard Deviation: $(diag(sqrt(stats_binomial.cov[:, :])))\")"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Calculate the cross section by sampling the same number of events from the efficiency and from the signal samples:  \n",
    "(Hint: proceed like in 4 a) for the samples of $\\lambda$_SB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resampled_SB = #=FILLME=#\n",
    "resampled_binomial = #=FILLME=#\n",
    "length(resampled_SB), length(resampled_binomial)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ_SB = valshape(prior_SB).(resampled_SB).params.λs\n",
    "ϵ   = #=FILLME=#\n",
    "σS  = #=FILLME=#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use a `StatsBase`histogram to visualize the cross section distribution.  \n",
    "From the plot, determine the 95% upper limit on the cross section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_σ = fit(Histogram, σS, nbins=200)\n",
    "plot(hist_σ, 1, seriestype=:smallest_intervals, xlabel=\"\\$\\\\sigma_s\\$\", ylabel=\"\\$P(\\\\sigma_s)\\$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the distribution of the signal cross section to the distribution from 4 a)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "hide_input": false,
  "kernelspec": {
   "display_name": "Julia 1.3.0-rc4",
   "language": "julia",
   "name": "julia-1.3"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}