Bootstrap confidence region.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## D. Heslop, J. L. Scealy, A. T. A. Wood, L. Tauxe and A. P. Roberts (2023). A bootstrap common mean direction test. Journal of Geophysical Research.\n",
    "### Example estimation of bootstrap confidence region for a mean direction (see Section 2 of manuscript) using the normal polarity directions of MacFadden et al. (1990).\n",
    "\n",
    "Please note that because this is a numerical procedure the results may differ slightly from those in the paper."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-W- cartopy is not installed\n",
      "    If you want to make maps, install using conda:\n",
      "    conda install cartopy\n",
      "-W- You haven't installed a module for plotting maps (cartopy or Basemap)\n",
      "    Recommended: install cartopy.  With conda:\n",
      "    conda install cartopy\n",
      "    For more information, see http://earthref.org/PmagPy/Cookbook#getting_python\n"
     ]
    }
   ],
   "source": [
    "import pmagpy.pmag as pmag\n",
    "import pmagpy.ipmag as ipmag\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```form_Mhat``` estimates $\\hat{M}$ according to equation (4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def form_Mhat(mhat):\n",
    "    \n",
    "    #input - mhat, estimated sample mean direction\n",
    "    #output - Mhat matrix according to equation 4\n",
    "    \n",
    "    b = np.matrix(mhat[0:2][:,np.newaxis])\n",
    "    c = mhat[2]\n",
    "    \n",
    "    if c==0:\n",
    "        A1 = np.eye(2)-(b*b.getH())\n",
    "    else:\n",
    "        A1 = c/np.abs(c)*np.eye(2)-c/(np.abs(c)+np.abs(c)**2)*(b*b.getH())\n",
    "        \n",
    "    A2 = -b\n",
    "    Mhat = np.hstack((A1,A2))\n",
    "    return Mhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```form_Ghat``` estimates $\\hat{G}$ according to equation (5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def form_Ghat(X,Mhat):\n",
    "    \n",
    "    #input - X, collection of directions (one per column)\n",
    "    #input - Mhat, Mhat matrix for mean direction\n",
    "    #output - Ghat matrix according to equation 5\n",
    "\n",
    "    n = np.shape(X)[1]\n",
    "    term1 = np.power(np.linalg.norm(np.sum(X,axis=1)/n),-2)/n\n",
    "    X = np.matrix(X)\n",
    "    Mhat_T = Mhat.getT()\n",
    "    Ghat = np.matrix(np.zeros((2,2)))\n",
    "    \n",
    "    for u in range(2):\n",
    "        for v in range(2):\n",
    "            for i in range(n):\n",
    "                Ghat[u,v] += Mhat_T[:,u].getT()*X[:,i]*X[:,i].getT()*Mhat_T[:,v]\n",
    "            Ghat[u,v] *= term1\n",
    "    \n",
    "    return Ghat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```find_T``` estimates $T$ according to equation (6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_T(m,n,Mhat,Ghat):\n",
    "    \n",
    "    #input - m, direction under consideration\n",
    "    #input - n, number of observations\n",
    "    #input - Mhat, Mhat matrix for mean direction\n",
    "    #input - Ghat matrix representing covariance\n",
    "    #output - T value estimated from Equation 6\n",
    "    \n",
    "    m = np.matrix(m[:,np.newaxis])\n",
    "    \n",
    "    return np.array(n*m.getT()*Mhat.getT()*np.linalg.inv(Ghat)*Mhat*m)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function ```find_CR``` calculates points along the boundary of the confidence region as described in Appendix A "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_CR(mhat,Mhat,Ghat,n,Tc):\n",
    "    \n",
    "    #input - mhat, mean direction of original data set\n",
    "    #input - Mhat, Mhat matrix for mean direction\n",
    "    #input - Ghat, matrix representing covariance\n",
    "    #input - n, number of observations\n",
    "    #input - Tc, critical T value on confidence region boundary\n",
    "    #output - mCI, closed confidence region boundary\n",
    "    \n",
    "    C = n*Mhat.getT()*np.linalg.inv(Ghat)*Mhat\n",
    "    [D,V] = np.linalg.eig(C)\n",
    "    \n",
    "    idx=np.flip(np.argsort(D))\n",
    "    D = D[idx]\n",
    "    V = V[:,idx]\n",
    "    \n",
    "    mCI = np.zeros((3,201))\n",
    "    y = np.matrix(np.zeros((3,1)))\n",
    "    for i in range(201):\n",
    "            theta = i*np.pi/100\n",
    "            y[0] = np.cos(theta)*np.sqrt(Tc)/np.sqrt(D[0])\n",
    "            y[1] = np.sin(theta)*np.sqrt(Tc)/np.sqrt(D[1])\n",
    "            y[2] = np.sqrt(1-y[0]**2-y[1]**2)\n",
    "            mCI[:,i] = np.ndarray.flatten(V*y)\n",
    "    \n",
    "    mCIbar = np.mean(mCI,axis=1)/np.linalg.norm(np.mean(mCI,axis=1))\n",
    "    if np.arctan2(np.linalg.norm(np.cross(mhat,mCIbar)),np.dot(mhat,mCIbar))>np.pi/2:\n",
    "        mCI *= -1\n",
    "        \n",
    "    return mCI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import the normal polarity directions from MacFadden et al. (1990) and transform to Cartesian coordinates (1 direction per column). \n",
    "\n",
    "#### Then estimate the mean direction via equation (3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = \"MacFadden_normal.csv\" #normal polarity directions of MacFadden et al. (1990)\n",
    "DI = pd.read_csv(file,header=None).to_numpy()\n",
    "X = np.transpose(pmag.dir2cart(DI))\n",
    "n = np.shape(X)[1] #number of directions\n",
    "mhat = np.mean(X,axis=1)\n",
    "mhat /= np.linalg.norm(mhat) #estimate mean direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main bootstrap sequence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Critical T value = 8.32\n"
     ]
    }
   ],
   "source": [
    "B = 10000 #number of bootstrap iterations\n",
    "T_b = np.zeros(B) #predefine output array for bootstrapped T values\n",
    "\n",
    "for i in range(B): #loop through boostrap iterations\n",
    "    idx = np.random.randint(0,n,n) #select observation indicies with replacement\n",
    "    X_b = X[:,idx] #form bootstrap sample\n",
    "    mhat_b = np.mean(X_b,axis=1) #mean direction of bootstrap sample\n",
    "    mhat_b /= np.linalg.norm(mhat_b)\n",
    "    Mhat_b = form_Mhat(mhat_b) #\\hat{M} for bootstrap sample\n",
    "    Ghat_b = form_Ghat(X_b,Mhat_b) #\\hat{G} for bootstrap sample\n",
    "    T_b[i] = find_T(mhat,n,Mhat_b,Ghat_b) #T for bootstrap sample\n",
    "    \n",
    "alpha = 0.05 #work with 95% confidence region\n",
    "Tc = np.quantile(T_b,1-alpha) #find 1-alpha quantile of T\n",
    "print(\"Critical T value = {:.2f}\".format(Tc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the 95% confidence region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mhat = form_Mhat(mhat)\n",
    "Ghat = form_Ghat(X,Mhat)\n",
    "Xc = find_CR(mhat,Mhat,Ghat,n,Tc) #calculate points along boundary of 95% confidence region\n",
    "DIc = pmag.cart2dir(np.transpose(Xc)) #convert to dec and inc for plotting "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the 95% confidence region (compare to manuscript Figure 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd7zcZZXH8c+XEEroHaUuRcAVpFioUqPoCiK6iC5NLOwqoigWdMXAKiiIKDbE3lY6AgrSVBAULFQRlBBpCoROICEk5Owf57mbyyXllvnNM+X7fr3mNblz5zdzkszMmef5Pc85igjMzMw6zSK1AzAzM5sXJygzM+tITlBmZtaRnKDMzKwjOUGZmVlHcoIyM7OO5ARlZmYdyQnKzMw6khOUmZl1JCcoMzPrSIvWDsCsV0haBBgHzA7XEDMbM/l9ZDaXpEWBNYB1gDWBVctlNWAVYNlyWaZcLwGMJ7/sjR/0UM8As4FZwHRgGvBEuX4MmFouD5Tre4C7gAciYk6Tf0ezbuEEZX1HkoBXkglmY2CTcv0vZFIaN+SQ2WQSeRB4nLmJZhqZfGaXyzuBlYFTyCQ0vlwmMDepLQOsQCa9lQENea6ZwN3AHcBtwK3AFODeiLitBX99s67hBGU9rSSj9YGXA1sCWwDbkEljwGNkIriDHMXcWa7vJUc4jw5nyk7SrwEiYqdhxjYOWAlYHVgLWJccua0LbEAmzSUHHfIgcA1wPXAdcG1E3D+c5zLrRk5Q1lPKFN1WwPbAduWyavn1M8DNwP3AvwFHAD8CprbinNFIE9QwHm8RMnF9Cng7cAawKbARcxc4TQGuLpcrgdt8/st6hRdJWNeTtCEwsVx2BpYrv5oCXEx+eF8L3BIRsyStQk7ZRUQ80MJQbmjhY1HORd0laQng7oh4C4CkpYCXkiPB7YDXAPuXw/4p6VLgMuAyj7Csm3kEZV2njJK2AfYslxeVX90FDHw4XxkR9y3gMe4HLoqItzcc7phJuhm4KyJeP5/fD0xj7kQm6V3JqUOA3wPnl8ufPbqybuIEZV2hJKWdgH2BvcgP4FnAL4ELyJHSHcP9AJZ0Cjm1d1QjAbdIST4/Ic83nTTMYxYBNgdeC+xBLgiBPLd2JnAacL2TlXU6JyjrWOXD+V3k+aKtyXNJ08jRwE+BSyLiiXoRdgdJLyD/Dd8IvJqc2r+dPB/3zYj4RcXwzObLCco6jqTVyXMqB5Mr2QI4i/zmf1FEzKgYXleTtCKwN7AfsGO5+Wrg28CZEfFkrdjMhnKpI+sISrtKOpdc3n088DDwQ3Kv0HERcU6nJydJk2rHsCAR8UhEfAv4cbnpy+R+rO8A90s6VdJm1QI0G8QjKKtK0jLAAcCh5GjpIfLD8jsR8dfyjf9+4CsR8cF6kQ6PpIiIoZtvO46kK8nKGC8uN21DbjTel9x7dSWZvM6LiFlVgrS+5wRlVZRpvPcD7yGrLPyR/EA8IyKeHnLfc4BtgTUjYna7Yx2JbkhQktYF/g58IiKOHfK7Fcmp1feSG4bvBr4AfCsinmproNb3nKCsrSRtQG6QPQhYjDy39AVyldo8X4yStiBr3l3TiSvPyrTep+bxq6MjYlJ7o1k4SePJvVPXRcQ/53OfcczdzLwD8Aj5BeIrEfFQu2K1/uYEZW0haT3gk+R03izg+8DnI+L2qoG1WDeMoEZK0rbAR8k9Z08CJwMnRsQjVQOznucEZY2StDbw32SpntnA14Hje7XCQS8mqAGSXgwcBbyFLJh7EnBSRDxeNTDrWV7FZ42QtIKkE8j9NgeSiWn9iPhgryan4ujaATQlIv4SEfsCm5HVOj4F3CHpMEmL1Y3OepFHUNZSpU7cYeS5ixXIqbyjIuKeqoFZy0naktwOsCswGTgO+G4nnie07uQEZS0j6bXAucDi5DfsIyLixrpRWZNKtY/dgc+TS9bvByZGxJ+rBmY9wVN8NmaS1pZ0NnAheRId4NNOTr0v0kXA98pNE4AbJJ0oadl6kVkvcIKyUZO0qKSPkJ1fXwt8nGy0dzfwpbJU2XqcpFXJhTAXAeuRZZMOB26T9O9llGU2Yk5QNiqSXkr2WPoccAmwSUQcFxGPAR8m+xUdXDFEa5//IUdOh0fEwxFxCFlB/T6yyeLZpWCt2Yj4HJSNiKTFyf1MHyU3b74XOHvwifHyjfmDwPe9qbP3SdoI2DYivjvk9kXJkdQxwEzyNeFFFDZsTlA2bJI2JVukb0auzvugN2s+l6RJnVg9oqbS8fhbwKuAnwHvbHEnY+tRTlC2UKUB3gfIZcSPAe+IiJ/Vjaoz9fJG3bEor6FDyWXpT5CvoQvqRmWdzuegbIFKUdeLgROBXwCbOjnZSEXEnIg4GdgK+CdwvqRTJC1ZOTTrYE5QNl+SdgFuALYDDgH2ioipdaPqPJImSQpJUX6OcplUObSOExG3kAsoTiBfU78rU4Bmz+MpPnuesjz8E2Qpm78B/+6Nl8PjKb7hk/RvwA/IqvbvjIjTK4dkHcYjKHsOScsDF5A15X4CvNzJyZoQET8HtgBuBk6TdFJZ+WcGOEHZIJI2AX4PTAT+C9g/Ip5c8FE2RM8Wi21CRNwN7Ei28PgA8AtJK9WNyjqFp/gMAEnHkHtWpgNvioirKodkFUlajdzvthZZwurUpvcvSXo7cArwEPC+iDinyeezzucEZUj6GLmEfAawRUT8tXJIVlGZ5r0FWAUYDzwFnBIRR7Thud8EnAkE8LqIuLjp57TO5Sm+PiZpnKQvkcnparIK+fGuodf39gCWJZMTwFLAYWUvU2MkrQB8htxr91fg55Le2eRzWmdzgupTkiYAZ5O9m04izwO8n2zrfVzF0Ky+8cDQlYiax20tI2k8WbdvPeCNwNZky5ZvSvqMC872JyeoPiRpOXLz7Z7AYaXL7bPAV4Gvlfv4A6F/XQjMAuaUn2cAp5fXSFMWA54GDomIKyLiCXIk902ySv6pHtn3H5+D6jOlNcLFwL8C/xERZw75vVzM00oB2C8CLyQriHwyIp5p+Dmf99orX5SOIdt5nEmuLJ3ZZBzWOZyg+oiktYFLyZVZe0fELyqH1HNcLLYZkj5Iltu6hHztPlU5JGsDJ6g+IWkd4FfAisDrvYy8Ga4k0RxJB5NTflcAezhJ9T4nqD4wJDntFhF/rBxSz3KCapak/yDLIzlJ9QEvkuhxTk7Nc7HY9omIHwMHkKtOfyZpqcohWYM8guphpc32b4CVcXJqC4+g2mPQSOpyciTlhRM9yCOoHlU2PV4MrA7s7uRkvaSMpN5F1o38sZeg9yYnqB5Upj1+DmxE9nC6pnJI/cTFYtskIr4DfAh4E/AN793rPZ7i6zFlR/75wKvJPk4uuGk9TdL/kPukPhsRR9aOx1rHvVd6SPkG+VVgd+DdTk7WJ44iC9t+TNLfI+LU2gFZa3gE1UMkHQkcCxwbEZ+oHY9Zu5RGhwMzB3tExEWVQ7IWcILqEaUq+WHA6cDbImLOQg6xHiRpcWAdYGpEPFY7nnaStAxZlf/FZBkvt5Dvcl4k0QMkbQm8p/y4LjndYX1G0suAfwJ/Au6X9N7KIbXb+sBK5OfayZL8PuhyTlBdrhR//SlwH/B2YFPg95I2qxqYtZSk90t6QNJDpf3EIkN+L3Ll5orA0szt7bVphXDbTtIbgKvIRof7A8sBZ5RFQ9alnKC6mKTFgLPIjbh7RcT3gB2AccDrK4ZmLVQ2pR4LrEqOED4ADO1uu3y5DPYs+YWlp5XkfCDZBfjlg/ZI7QR8oWJoNkZOUN3teDIhvSMirgMo15vjpoNVNFTe6C3AhEE/TwD2HXKfx4Gh1RQWAe5oIJ6OUlp0HADsFBH3ldt+SCanQyXtVzM+Gz0nqC4laW+yA+7JEfGTwb+LiIfc06maTzXwmI+Qo6EBUW6be0Muinkz8BSZrGYAX4mIaxuIp+NExJMRMWPIzR8lS319Q9ImFcKyMfIqvi4kaT3gOuBvwPZNN5Kz4WuiFp+k9cmFDxPIL5VPAztGxJ/mcd/VyWaU/4yIW1sZRzeStAZwPTAVeKWrn3cXJ6guU847/ZZcsbRFRNxZNyIr03rzGjkd3armhZLWAvYjzy+eGRF/bcXj9gNJrya7An8vIg6uHY8NnxNUl5F0HPAxsqvoubXj6XaSNga2Bh4ALh7r/jFXM+9Mkj4NfIIs/3VW7XhseJyguoikVwG/Br4dEe+qHE7Xk/RG4EfkOZ0AriSrEIw6STlBdaay3Py3wHrAphHxz8oh2TA4QXUJScsBNwKzgc0j4snKIXW1sjT5cWCZQTc/Cbw1In62gOOWA2ZFxPT5/H5Sq6b1rLUkbUSej/oN8FpXW+l8XsXXPb4ArAns5+TUEouSG1oHE/DCed1Z0tKSfgU8CDwu6evzau/g5NS5ynm7D5L1+v6zcjg2DE5QXUDSROBg4AT3dmqNiJgF/JnnLt9eBJjfv++XgW2A8WRyOwB4Z5MxjobShpK2kLRE7Xg60DeAS4HPSVq7djC2YJ7i63CSliY/SJ8mp/aerhxSzygfUBcCm5CbXN8xdE/ZoPtOJldODvaTiHhbs1H+fwPKg8l9ThsDy5KJcmAEF2SinUEm2SWBWcDDwA5e6flcktYl31NXAv/mPYOdyyOozvcZYG3yw9PJqYUi4u6IeAmwFLDU/JJTcRcw+JzFTGByU7GVkdCxkqaR58ZOBl5Fljtaglxuvki5jAMWI+vPLUOO8JYkp4QnSzpf0uZNxdptSsI+EngtuXTfOpRHUB2sVKf+PfC1iDi0djz9TNKGwO+YO3K5B9g6IqY18FwfI7+YDP4C+SxwN7kS7UpyBHB/uX0pMhl9gkxi8zOdbGj58YiY3eq4u0kptvsb4EXARhHxyEIOsQqcoDpUeQP9juzts1FEPF45pL4naUUyAcwEftXEiFbSbcBGg256FphCVkF4dCHHHkyOtJYadOw95IhqzSGPeSpwaD+vZCsV/68DTo2I9yzs/tZ+nuLrXAcDrwCOcHLqDBHxSET8tHRrXVHStZJmlTYYrxnr40v6G3OT08A3x3Hkl5Tjh/EQ3wMuIs9FTSNbsOwUEWuR561+QG5TGAf8FzBdUt9WvY+Im4CvAP8paava8djzeQTVgUrttVuAm8hvzv5P6jCSbiQ7ty5abpoObAbsP5ql5pIOAr5bfnwUWGHIXX4fEa8cxuMI2IBcQn/r0FFeGZmfQq5AHFhk8dOIeONIY+4FZV/bFHJU+YKIeHYhh1gbeQTVmY4nG85tDhxXWllbh5A0gSzIuuigm58FtmX01cy/Wq5nAd8iR0EDngb+MJwHiXR7RFw/rynIiJgTEe8muy7fU27eS9Id89rX1ctKfcOvk00eVwHeVzciG8oJqsNI2gDYgyzB8xOyZcDt7pDbUZ4mp8oGE7mse7SWLNf/C0wC/kiOyqaTFUQ+NobHfp6IeDgi1iZfZ5AlgP7eyufoZJLeBPwVeCPwaeAG4HBJSy7wQGsrJ6jOcxzwDPDhiDgQeCVZf8/VqyuStJikN0p6O3lO6DAyecwk/7+WJluuIynKZdJInqJcP1HKKO0IvBTYAth2ftVDSoWLkyX9TtI3Ja0r6XRJd0q6vLRmma+I2B/4fPlxHUnnjyDmUSlL6Fcq02u1/An4KbBJRHySrDCxNvn/ah3C56A6iKRtgauBoyLif2rHY6lUZLiKuQsYBLyOTEzbAP8EzoqIZ0dbLFbS0+S07mMRMfT80/yOWaTEtQW5N2omubhiEXJf1LPkqO5FC1toI+kyYNfy4wYR0Ugn3rLx/AJyOlTA6cBBnXDupyTnHcm//4O14zGPoDrNp8m2D1+oHYg9x0FktYmly2Up4LsRcU1EnBQRp7fgA/akcr28pL2Gecx65ChroKTR4uXPi5Wfx5Wftx7GY01kbtmny4f5/KNxIhnPYuSesr0poxZJ60l6v6T/lLRSgzHMz0fJ/98PV3humwcnqA4haWdgZ+A4d/3sOKsz9xzRgFXmc9+jR/MEEXEkOQICOHeYlR+GM/0h8pzZwp4/yHOekFNdTdmBuQkVskvwDmWZ943AZ8kvaLdIWq3BOJ6ndCD+X+DQ0pnYKnOC6gBl9dQx5FTRNyqHY893Fc9dVfcMORX7PBExSdJSo1wRtwFzk871khbW82sKeS5lILaZ5BL1gS84M8hzl/OMdR4+Uq6bXM13B89dYDITuB34EjkyXYL8MrAiOaJpt2PI0V2N57YhnKA6w27A9sBnXG+v80TEZcDHycQ0h1zy/R9D7ydpI0lTyD5TT0m6QtIPhrsJNCLuJfdWDUy1nSrplvktJiijnt3JJepXksvT1yO3KVwDnAfsOoKyRkOL4TbhveR5sSfIzcRTyLJOq/DcxDgeeEEb4nmOiLid3ND8X5Lm2XrF2seLJDpA6TO0IbB+RMxc2P2tjjIqGh8Rz8znd3eRJYWGjkCmAztGxB+H+TwTyJHPQHmiAC4G9h3GYofXAOeUH+eUx9luOK8rSQ+TI5c5ETFuOLGOhqRlyS9ks4ErImKmpBOA95BTfpD/ZgdHxOlNxbGA+NYD/gacFBE+H1WRR1CVSXolsBPwBSenzlY2wT4vORUrkZXG5zU9NgE4YgRPtRhzqxtQHnN34DFJt0vaawFTiN8uzzeBPOG/McOo2C3pSjI5QZZMWtB9d5X0JUlHSVp5YY89VEQ8EREXRsQlg17zHwdOI8+XTQOOqZGcSnxTgDPIEkjDWlFpzfAIqjJJ55CLI9ZuojK2tYek8eQH6+Lzuct5ETGs1XmSziWXsQ+sxnuGTFLjB91tDvkt/wfANyPioXLsdJ67oGM2uW3huHk8j4BDyZV1A4/9UETMbwEIkg4gqy9MIKtePARsGhFj2aTcccoileuBT0TEsbXj6VdOUBVJ2gi4lTz39Mna8djYSDoQ+Bo5MzF4pdp0YO+IuHiYj/MQOSIb7PPk8u+TyBYRQ2c/niUXHMwpzz1QhmkGWXj4VrKd/RZkRfaXkSOmwSOxm8immPP9UJB0H7mqccDTwJER8cXh/N26iaSLgC2BdXxuuI5FF34Xa9D7yG+hX64diI1dRHxf0nXkh9omwGvIEcxnhpuciqk8N0HNAO6JiF9I+i1ZSHg1njuiGsfc8zeDLcnc5ePzcx+wV0T8fhixDV1uvyhz23v0mhPILwX7spBpT2uGR1CVlBPF/wDOKSWNzACQtDVwafkxyPNR20TEjDLF9lXy/NKAmWR32MOBNcjVb8uRCWzwYocgvxA9AfyFrJ7+45FsMpb0beCtzE1U08lSTDeO5O/YDcoU6M3kv+/L3FWg/TyCqucg8kPGo6ceImnSaNptDBYR10j6V7LszpPAhYMWEyzO86f3xpGr4X41JJblgH8hR1+tOkf0HvIDey9yOf37ejE5QS6KkfQV8pzbNmQ3Y2sjj6AqKDXUbgUeiYhtasdjrTPaWnwjePy1yCm+pcnzRzOAn0XEPkPu9zpyJdqz5GKL/4yI7zcVV68qtQPvBS6KiLfWjqffeJl5Ha8iT3R/rXYg1l0i4h5yD9HVwGRyc+4Bg+9TPlTPIM8NLUsumvh6SW42AqWK/PeBvSvVB+xrTlB1HEyeBzi7diA2dpImDbTYKD+Ppt3GsEXETRGxQ0RsGBGHzWOF2Zo8v07fM+RmcBu5b5Oj0OdVD7FmeYqvzcp5gfuA70fEf9WOx1qr6Sm+YcawDHA/z13VN4PsfXRXnai6m6Q/kElqgcvwrbU8gmq/z5EroBpvDGf9qWz43o9cYfc4mZw+6OQ0JhcCm5FLzq1NPIJqs0EbHeeQeyx+SDa7m7HAA60rtGIVX6tIWpWskH5nRPxzHr9flKxgvhtwJ/DxiLi/rUF2sFLGaT9gf3JvG8AvI2LX+R9lreQE1UblJPXdZDWAJ8kX/xrACyLikZqxWf+R9H3gzcwtW/QgOQ34RNXAOoSkHYFfky1NfkhuvN4cWKsTOgD3AyeoNpJ0BLk7fcOImFw2Am5cGqWZtY2kxcm+UYM38k4jK4ifVSeqzlK2g2w08P6U9BayoO3OEfHrmrH1C5+Daq99gT9ExGT4/+rYTk7WSaou8OgkETFnyPvzZ2RS93moNnGCahNJ/wJsRe5PMauqVKY4i1xIAVkzcDpzSyzZEBHxFHAB8CZJjfXLsrmcoNpnz3J9btUozObaH/gC2X33LLLe3GN1Q+p4PwVWBrauHUg/8DmoNpH0S2DViHhJ7VjMOlnpKPxKclR3TUTMqhzS/yv7GB8EvhgRH6kdT6/zCKoNSlfOV+G9T9YgSUtLerukwyRt3Obn3k3SnZKekHSBpOVH+TgvBG4jRyo/B/5YSjd1hIh4nFzZt+dC7mot4ATVHruTq6WcoKwRpXrE9WR1/M8Bf5K0c5ueeyPgPGAdYBlgIjllOBonk+1Cli2PtRHw3y0Is5XOAzaS5NJRDXOCao+JwKPAH2oHYs1qqv7eMBwCrEUWiF2C3Nv0jTY99y5Dfl4c2Kks0x6pjXluG6DFgRePNrCGDDSfnFg1ij7gBNWwstdpInC5N/f1hU9Vet5VyVpxg63Qpud+jKyMMtjTETH0tuG4huw3NWA6ndeH6Q6y8oYTVMOcoJr3IrK6tJfvWpMuIWvuDXi63NYO55Jdf6eTiWo62d13NA4npypnkInqEuDEFsTYMqVY7GXAzqVclDXECap5u5Xry6pGYY1pd7uNeYmIy4APkW1cniGnoQ5p03M/TS67/hBwNPDaiPjmKB9rGrAt+cVu3Yh4Yyet4hvkUmA54GW1A+llXmbeMEmnA1tHxDq1Y7HmdUK7DWteKcT7APCRiDihdjy9yiOoBpXzTwPdT82sR0TEVLKj8fa1Y+llTlDNWgd4IU5Q/eTo2gFY21wFbFu+iFoDnKCatV25doLqE53SC8ra4mqy7NGLagfSq5ygmvVKsvrxzbUDMbOW+125dl2+hjhBNWtL4AbvfzLrSbeRy+G3qB1Ir3KCakjZRb85cF3tWMys9coXzxtwgmqME1RztiLLztxYOxAza8xNwFbesNsMJ6jm7FeuT5F0s6T/lbR71YjMbMwkvVvSuZImA+8mv4juWDmsnuQE1ZzZQAAnAXeRK/rWrxqRmbXCK4BNyJJM3yq3rVQvnN7lYWlz1gKmDG5q5v0SZj3hXaUe30Cvt3cB61aNqEd5BNWcjYFbB98Qritl1vUGv48j4lGy5FFbG0T2CyeoBkgaR27eu3Vh97XeUrEflNVzKznlZy3mBNWMNchGa5NrB2JtV6sflNUzGZ9fboQTVDMGKpffWTMIM2uLO4FVJE2oHUivcYJqxkCCuqtqFNYWndAPyqq6s1yvXTOIXuQE1Yx1y/XdNYOw9oiISRGhgT5QA3924di+MfBFdN2aQfQiJ6hmrAU8GBEzFnpPM+t2A19E16oaRQ9ygmrGqsD9tYOwKtwPqv88UK5XrRpFD3KCasZqwNTaQVj7eVqv/0TETOBx8n1vLeQE1YxVcYIy6ydT8Qiq5ZygmuEEZdZfpuIRVMs5QbVYqSKxDPBo7VjMrG0eBZarHUSvcYJqvaXL9bSqUZhZO00jv5haCzlBtd6y5doJyqx/TGPue99axAmq9Qa+RT1RNQoza6cn8Aiq5ZygWm+gHtf0qlGYWTtNB5Z0z7fWcoJqvfHlelbVKMysnQbe724C20JOUK038AIddoIaaVHR0RYh7Zbjah7bzc891uNrHNvO4xo+xgmqAXKT19aStAtwObBTRFwxzGNioNBoE/fvtuNqHtvNzz3W42sc287jmjxG0oeAzwPLRYTPP7eIs33rHVquT5I08EK9ISI+UCsgM2s9SV8ENi8/rlGu/ZnaQp7ia72BIekCv3WNtIfQaHsOdctxNY/t5ueuGXs3vLbadQxz3+/PLujvYCPjKb4Wk/Qq4Apg14j45TCP8RRfhxzbzc891uM9xTf6YyR9BPgcsFREeAVvi3gE1Xqzy/X4Bd7LzHrJwPt99gLvZSPiBNV6A6t5RpKgRtpDaLQ9h7rluJrHdvNzj/X4Gse287gmjxk49+QE1UKe4msxSZsCNwH7RMSZteMxs+ZJ+ixweEQsXjuWXuIRVOsNrNxz2ROz/rEsLm/Wck5QrTdQJNYJyqx/LIMLRLecE1TrDbxIXdnYrH94BNUAJ6gWi4hZZOFINy8z6x/L4QTVck5QzXgQWKV2EGbWNquQbd+thZygmjEVWK12EGbWNqvhBNVyTlDNmAqsWjsIM2uepEWBlXCCajknqGY8gEdQfWm49fOspwxM5ztBtZgTVDP+AawuyeWO+s+nagdgbbdmub63ahQ9yAmqGXeR/7ZrLuyOZtb11i3Xd9UMohc5QTXjznK9Ts0grD3G2i7Dut665doJqsXcXKsZAy9UJ6g+EBGTgEkw9nYZ1pXWAR6PiMdqB9JrPIJqxt1k47INagdiZo1bH/h77SB6kRNUAyLiGWAKsHHtWKztxtouw7rPxsCttYPoRU5QzbkNJ6i+U6b7rE9ImkBO8d1WO5Ze5HNQzbkN2F3SPsBLgE2BH0XE2XXDMrOxkHQk+X6+GXgSEB5BNcIJqjkrkV11TwfmALfjCudmvWB5YDvgrYNue7ZSLD3NHXUbIul1wM+BjwEnR8SMyiGZWQtJWh74KrAPsFxETK8cUs/xOajmXAbMAlZwcjLrPWVZ+WrADU5OzXCCakhZyfdnYMvasZhZ60kS+f6+vnYsvcoJqlnXAVuWF7KZ9ZZ1gBXI97k1wAmqWX8gF0usXzsQM2u5rcv1H6pG0cOcoJp1dbnermoUZtaE7YCngBtrB9KrnKCa9RfgMZyg+oYLxPaV7YFrI2J27UB6lRNUgyJiDvBbnKD6iftB9QFJywCbAVfVjqWXOUE17yrgxZLcAt6sd2xHfn46QTXICap5l5XrXapGYY1xP6i+NBGYiRNUo1xJomGSxgEPAudGxDtqx2PNcj+o/iDpRuChiNi1diy9zCOohkXEs8AvgYneD2XW/SStRp5/urR2LL3OCao9LgXWwu03+oH7QfW+V5fryxZ4LxszJ6j2uLBc71E1Cmuc+0H1hT2A+3AFicY5QbVBRNxDvpj3rB2LmY2epMWB3YELyjYSa5ATVPucD2zr5eZmXW0nYBny/WwNc4Jqn/PIzpue5jPrXm8AprXgeFIAAA5tSURBVAOX1w6kH3iZeZuUFXyTgckR8Zra8ZjZyEhaFPgHcEVE7FM7nn7gEVSbRH4TOA3YdfA0X5nTNrMONOT9uTOwKvCTSuH0HSeo9joNGAfsI2kvSWcDd0taonJcZjaEpE2BBySdKmkHYF9gGnBR3cj6x6K1A+gzfwamAieR//ZTgf8FJgBPV4zLzJ5vFrkY4m3Au4A5wJUR4fdqm3gE1UZlmu9PZHJ6B7BGRBweEY/UjcxaxfX3ekdE3BYRBwCrA18mPy8vqRtVf3GCar93k9/E1nMfmZ7kdhs9JiKeBF4E3AMcXzmcvuIE1WYRcS9wMXBQKSRrZh1M0lpkeaPvldqa1iZOUHV8G1iDuTW9rIu53UbPO4jcw/jdynH0He+DqkDSYsC9wDUR4fJHPcTtNnpL2fs0BfhbROxWO55+4xFUBRHxDHAq8HpJ69WOx8zm6w1kJ4Kv1A6kHzlB1XMKuVjiPbUDsZZyu43e8j7gLuCC2oH0I0/xVSTpDLJ19JoR8VTteMxsrrJR9ybgoxHh1XsVeARV18nA8sDbawdiZs/zIWAGuajJKvAIqqJSQPZq4AXAht4XZdYZytLyKcDXI+Kw2vH0K4+gKiqVJT4LrAu4OrJZ5zicXFp+Yu1A+plHUJVJWgS4GZgNbB7+DzGrStKKwN3AuRGxf+14+plHUJWVttGfAzbDLeHNOsHhwFLk+9Iq8giqA5TNgH8hT8huUZKWmbWZpJWBvwMXRsRbasfT7zyC6gBlccTR5Chq78rhmPWzI8jRk/ezdQCPoDpEKRx7MxDAZi5K2Z0kTYqISbXjsJGTtBq5cu/ciNivdjzmEVTHKAnpU8CLgQMqh2Oj53Yb3esoYHHgmNqBWHKC6ixnAdcCn5a0FICklSR9uOyZMrMWkLS8pA9JWqL8vAlwCHBKRPytbnQ2wAmqg5Ql5h8CXggcIekw4HbgOGCrmrHZ/LndRld6M/B54BZJe5Er9qbjc08dxeegOpCky4GdyY2ClwGHR8Sf60Zlw+F2G91D0kTgi+S0OsCxEfGJiiHZEE5QHUjSrsClwBXALt682z2coLpLmeKbAqwErB4Rj1YOyQbxFF8HiojLgc8AOwE71I3GRshTRN3lPWQtzH2cnDqPR1AdStIEcvPuNGDLiJhVOSSzniLphcBtwG+A13umovN4BNWhImI68AHgJcD7K4dj1otOBBYDDnNy6kxOUJ3tPLKT5zGS1q8djFmvkLQHsC+5MOKO2vHYvHmKr8NJWoOc6vsTsKu/6ZmNjaTlgFuAR4CXRcQzlUOy+fAIqsNFxD/I+mA7A++sHI5ZLzieXBjxDienzuYRVBcoVSQuB15G1um7s25EZt1J0quBi4HPR8SHa8djC+YE1SUkrQvcBNwA7OxismYjI2klsiDzo+TU3ozKIdlCeIqvS5RR06Hkvqgj6kZj8+PyRp2pzEJ8A1gZ2M/JqTt4BNVFypvsDOANwNYRcV3lkGwIV5LoTJIOBL4HfCwi3Cm3SzhBdZkyTXEj2X13q4h4onJINogTVOcplcr/WC67eHq8e3iKr8tExMPk/o1/Ab7lNhz1uZp55ypta84EngLe6uTUXTyC6lKSPgp8Fjg0Ir5aOx5LHkF1FknfBQ4EXhMRl9aOx0bGCapLSVqErDKxG7BjRFwz+HcRMadacH3MCaqO8n6IwRvZJb0T+CbwPxFxVLXgbNQ8xdelSgI6APgHcE4pfImkVYArJP17zfj6mKuZ13Es8B1JiwNI2hb4GnAJ/j/pWk5QXaycj3oDsCyZpLYEfk9u6PXQuIKImFQ7hn5TzsM+DRwEXCZpM+Ac4G5gX5936l6e4usBkt4EnAXMAh4G9oyIP9SNyqy9JO0LfJfsRP0s8IqIuKVuVDYWHkH1gIg4G7gSGA9838nJ+tQZwG+BxYGTnZy6nxNU79gZOA34qKS31Q7GrIITgF3IzbhH1g7Gxs5TfD2knCC+GNgGeHVEXFE5JLO2kHQo8OVyeb/b0vQGJ6geI2kFcprjBcBOEXFD5ZDMGiVpH3L24AJgby+K6B1OUD1I0trA1WQ76+0j4vbKIZk1QtLuwPnAteRm3OmVQ7IW8jmoHhQRdwMTyf/fyyStWTmkvuHyRu0jaXtyOfktwB5OTr3HI6geVvZF/Qq4j+whdV/lkHqeK0m0h6SXA5cCDwA7RMTUyiFZAzyC6mGlHce/AWsCv5L0gsohmY3ZoOT0CDDRyal3OUH1uIi4CtgdJ6nGuJp5+wxJTjuX6WzrUZ7i6xNlvv4XZO2+3SLinsoh9SRP8TWn1Ne7kLnJ6a7KIVnDPILqE4NGUqsDV0l6UeWQzIZN0qvJkdOD5PYJJ6c+4ATVR0qS2glYEviNpM3rRtSTXDm7xSS9GfgZcDu5bcLTen3CCarPRMT1wA7ATODXknYe/HtJh5SimzYKrmY+OpIWkXS0pB2G3H4IcDpZpX+niHigSoBWhRNUH4qIvwLbk+ejLpa0P4CkicBXgbe5lby12dLAvsC5ktYrCetzwCnAReQm3MeqRmht50USfUzS8uRGx53JGmYHkD10touIaTVjs/4jaUOyIsQDwG3AXsDXgcMiYnbN2KwOJ6g+J2kx4HvAW8mmby+NiL9VDcr6VultdibZ0+nDwIku/Nq/nKBsoCPppcCuwB+BN0bEvXWjsn4j6RXAucDKwJcj4ojKIVllPgdlRNqNbB+/EfBHSdtVDsv6iKQDyaabzwAvd3IycIKyQSLifGBrYBq5wu9wL5awJklaQtIp5DTz1WRyuqluVNYpnKDsOSLiL8DLgZ8DXyBXVa1QN6ru4fJGw1cWRVwDHAIcT67Ue6huVNZJfA7K5qmMnD5AfnDcC+wbEdfWjarzudTR8JS9dqcCs4ADI+JnlUOyDuQRlM1TOS91ErlfSsDVkj4ladHKoVkXk7ScpB8CPyH7OG3h5GTz4wRlC1RGTS8lP1AmkXX8NqgaVIdxNfPhkbQjcBO5pWES2cfJZYtsvjzFZ8MmaR9yZ//iwJHAVyPi2bpRdRZP8T2fpKWATwPvByYD+3u62IbDIygbtog4A9gUuAL4EllwdpN53VfSzpLGtTM+q0PSUqUVxrx+txvwZ/J85tfJKT0nJxsWJygbkYj4B9mld39yz9QNkj4pafGB+5QPq18C760TZVX9WM3868CFklYZuEHSypK+Q24AfwZ4VUS8NyKeqhWkdR9P8dmoSVqVrOG3Dzl18z7gErKe2guBjSLiyXoRWjuUUfTNwLfILyXvBI4FlgNOAI6OiKfrRWjdygnKxqw0k/sy8CLgT8BW5HmGH1UNzNpG0hfJLyi3Av8K/Bo4NCJuqRmXdTdP8dmYRcQlwGbAJ4AtgQBeImm5qoFZW0haB1iD/DxZH3gbsIuTk42VR1DWUqXg55FkXb+HgWOAb0TEM1UDs5Yr7VqOJFfnBXAa8N/lPKXZmDlBWSMkbUmef9gF+DuZqH7kvj7dT9IywGHAEeR5ph+QickV8K2lPMVnjYiI64DdgNcBjwDfBf4i6T96efl5L2/OLcvJP0J+4fg08Btgy4g4yMnJmuARlDWu1PXbkxxFbUau+DsB+EGvre7qxY26klYCDiUXQawE/AI4KiL+UDUw63keQVnjSl2/84AtgDcBjwHfAO6U9LFyLmO+JL3QCy5aS9KaZapuQfdZp6zOu5ssTfRbYNuIeK2Tk7WDE5S1TUTMiYhzgFeQ3XtvBI4D/iHpFEmbzufQY4DJksa3KdQR6bZafKXNxd3karuhv5OkXSSdC0wh9zWdBbwkIvaMiN+1N1rrZ57is6okbU5OHb0NWILcP3MK8NOImClpCeB+4LyIOLBaoMPUDVN8Zcr1z8BjEbFduW0F8v/gPcCLyRWY3wS+FhH31IrV+psTlHWEcp7jYPIDcl1yYcWPgXvInlS7RcTl1QIcpm5IUACSPgp8FjgAeC2wN1kE+Dpy0/VpvXZ+0LqPe/tYR4iIh4ETJJ1ITv8dTHZaXYxsavcKSVMi4u8VwxyOjq7FV0ZPWwDrlJt+ADxKjpa+ExHX14rNbCiPoKxjSVoR+BrZgn69cvO1wBnkFOCUWrF1k0FJaS/gLWRJqtnAXcA55Io8j5as4zhBWVeQtC5ZlPatwObl5luA84ELgN+PpDeVpPcDj0TED1sbaWuV5HIk8NeIOHsExy0J7Egu798DWBOYA/yKrPhwbhm1mnUsJyjrOpLWIz949wReBYwjl67/imzvcBkwORbw4pb0d+DaiNi3+YjHRtJk4PqI+PcF3GcRMnFPLJftyXNK04GLyUT+84h4sPmIzVrDCcq6Wll99hqyasVEYO3yq/uAq8vlKuDGiJhVjlkGeIIsz/OZFsbyRYCI+ECrHrM87rnAxhGxyaDbJpDL9bcjk9E2ZNkhyBV6A4n6l56+s27lRRLW1SLiUXLK6rQyHbYBmay2Jz+831zuOlPSTcD15GgLsqJFK22+8LuMTPk7TQb2LGWGNiXPJ23C3H2MtwCnk6WHLo+I+1odh1kNHkFZT5O0BpmoXk62AtkCWGHQXf4B3Eb2MZpCLhy4s1w/sqBpwnk8168BImKnEca4CLAaubx+nXLZENiYTEQrDon3+nK5BvhdSdJmPccJyvpKGZG8jNz38wTZtn4TMhksO+TuM4Cp5fIA8CDweDluWrmeQa6Im0WWAwrgKHJ2Yny5TCiPvUy5LE8mpFUHXYZWyXiQTJq3knvBAH7ilYvWT5ygzPj/xLUCzx3FrMncBLIasApzE81ozCET2+Nkwps66PoeBo3eIuLJUT6HWc9wgjIboTIltxSZrJYgRz8DI6ZxzB1RzS6Xp8jENH0kU4Zm/c4JyszMOpKrmZuZWUdygjIzs47kBGVmZh3JCcrMzDqSE5SZmXUkJygzM+tITlBmZtaRnKDMzKwjOUGZmVlHcoIyM7OO5ARlZmYd6f8AGSDEJ/Hl1GcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ipmag.plot_net(1)\n",
    "ipmag.plot_di(dec=DI[:,0],inc=DI[:,1],color='black',edge='black')\n",
    "ipmag.plot_di(di_block=DIc,color='black',edge='black',marker='.',markersize=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Test if any known point is inside 95% confidence region\n",
    "If ```T0 < = Tc``` point is inside confidence region\n",
    "\n",
    "If ```T0 > Tc``` point is outside confidence region\n",
    "\n",
    "As an example case, choose a point on the 95% confidence region boundary and check that it returns ```Tc``` as calculated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New point T value = 8.32\n"
     ]
    }
   ],
   "source": [
    "T0 = find_T(Xc[:,0],n,Mhat,Ghat)\n",
    "print(\"New point T value = {:.2f}\".format(np.squeeze(T0)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}