{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "13fae235",
   "metadata": {},
   "source": [
    "### Calculate error statistics from array of weights"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c4fd2af",
   "metadata": {},
   "source": [
    "Import packages. This only requires jax to work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e8acfe1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "\n",
    "# load in population-error package\n",
    "from population_error import error_statistics_from_weights"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e05afcf",
   "metadata": {},
   "source": [
    "For each found VT injection $\\theta_j$ and each sample from the hyperposterior $\\Lambda_n$, `vt_weights` array should be $p(\\theta_j | \\Lambda_n) / p(\\theta_j | {\\rm draw})$ weights. The shape is `vt_weights.shape = (Nsamp, Nfound)`\n",
    "\n",
    "Similarly, for each posterior sample, for the $i^{\\rm th}$ event, the $j^{\\rm th}$ posterior sample $\\theta_{ij}$, `event_weights` should be $p(\\theta_{ij} | \\Lambda_n) / \\pi(\\theta_{ij}|{\\rm PE})$ are the weights. The shape is `event_weights.shape = (Nsamp, Nobs, NPE)`.\n",
    "\n",
    "Finally, `total_samples` is the total number of injections, $N_{\\rm inj}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "092790dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# just random arrays for simplicity. Often this approach runs into memory errors\n",
    "vt_weights = jnp.array(np.random.uniform(size=(500, 100_000)))\n",
    "event_weights = jnp.array(np.random.uniform(size=(500, 100, 1000)))\n",
    "total_samples = 5e6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f57db5e",
   "metadata": {},
   "source": [
    "Compute error statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ac52071e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Running for 500 iterations: 100%|██████████| 500/500 [00:10<00:00, 45.96it/s]\n",
      "Running for 500 iterations: 100%|██████████| 500/500 [00:10<00:00, 45.94it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'error_statistic': 0.14175374994676887, 'precision_statistic': 0.14175374564949028, 'accuracy_statistic': 4.297278581024669e-09}\n"
     ]
    }
   ],
   "source": [
    "statistics = error_statistics_from_weights(vt_weights, event_weights, total_samples, include_likelihood_correction=True)\n",
    "print(statistics)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gwjax311",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}