PX_prior_testing_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from cycler import cycler\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.stats as sps\n",
    "import json, pickle, os, corner, glob, sys\n",
    "from collections import defaultdict\n",
    "import enterprise.signals.utils as utils\n",
    "from enterprise.pulsar import Pulsar\n",
    "from collections import defaultdict\n",
    "import pandas as pd\n",
    "\n",
    "import enterprise\n",
    "from enterprise.pulsar import Pulsar\n",
    "from enterprise.signals import utils\n",
    "from enterprise.signals import parameter\n",
    "from enterprise.signals import white_signals\n",
    "from enterprise.signals import signal_base\n",
    "from enterprise.signals import selections\n",
    "\n",
    "from PTMCMCSampler.PTMCMCSampler import PTSampler as ptmcmc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyne2001\n",
    "import scipy\n",
    "from scipy.stats import rv_continuous\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_cycle_wong = ['#000000','#E69F00','#009E73','#56B4E9','#0072B2','#F0E442','#D55E00','#CC79A7']\n",
    "mpl.rcParams['axes.prop_cycle'] = cycler(color=color_cycle_wong)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "current_path = os.getcwd()\n",
    "splt_path = current_path.split(\"/\")\n",
    "#top_path_idx = splt_path.index(\"akaiser\")\n",
    "top_path_idx = splt_path.index(\"nanograv\")\n",
    "top_dir = \"/\".join(splt_path[0 : top_path_idx + 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "e_e_path = top_dir + \"/enterprise_extensions/\"\n",
    "sys.path.insert(0, e_e_path)\n",
    "sys.path.insert(0, top_dir + \"/enterprise_timing\")\n",
    "import nltm_plot_utils_v2 as nltm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/andrew/Documents/Research/nanograv/enterprise_timing/J1640+2224/chains/\n"
     ]
    }
   ],
   "source": [
    "psr_name = 'J1640+2224'#'J1640+2224'#'J2317+1439' 'J1909-3744'\n",
    "psr_dir = top_dir + '/enterprise_timing/{}/chains/'.format(psr_name)\n",
    "psrlist = [psr_name]\n",
    "datareleases = ['11yr','12p5yr']\n",
    "ephem = 'DE405'\n",
    "\n",
    "#psrs = nltm.get_psrs(psrlist,datareleases)\n",
    "#pardict = nltm.get_pardict(psrs,datareleases)\n",
    "print(psr_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = top_dir + \"/5yr\"\n",
    "parfiles = sorted(glob.glob(datadir + \"/par/*_nltm.par\"))\n",
    "timfiles = sorted(glob.glob(datadir + \"/tim/*.tim\"))\n",
    "parfile = [pfile for pfile in parfiles if psr_name in pfile][0]\n",
    "timfile = [tfile for tfile in timfiles if psr_name in tfile][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter\n",
    "is_psr = False\n",
    "if psr_name in parfile:\n",
    "    psr = Pulsar(parfile, timfile, ephem=ephem, clk=None, drop_t2pulsar=False)\n",
    "    is_psr = True\n",
    "\n",
    "if not is_psr:\n",
    "    raise ValueError(\n",
    "        \"does not exist in datarelease.\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('RAJ', 'DECJ')\n",
      "[4.36454075 0.39099711]\n"
     ]
    }
   ],
   "source": [
    "print(psr.t2pulsar.pars()[0:2])\n",
    "print(psr.t2pulsar.vals()[0:2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def px_prior(px,dist,dist_err):\n",
    "    y = 1/(np.sqrt(2*np.pi)*dist_err*px**2)*np.exp(-(px**(-1)-dist)**2/(2*dist_err**2))\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = SkyCoord(ra='16:40:16.7435004', dec='22:24:08.94331', unit=(u.hourangle, u.deg))\n",
    "c_gal = c.transform_to('galactic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41.05104415279213"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_gal.l.degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.16"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dm_dist = pyne2001.get_dist(c_gal.l.degree,c_gal.b.degree,18.426)[0]\n",
    "dm_dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NE2001DMDist_gen(rv_continuous):\n",
    "    \"Gaussian distribution\"\n",
    "    def _pdf(self,x, mu, sigma):\n",
    "        return (1/x**2)*scipy.stats.norm.pdf(x**-1,loc=mu,scale=sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_diste = NE2001DMDist_gen(shapes='mu,sigma')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.991870553452739e-07"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_diste.pdf(10,mu=1.0,sigma=.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8252921321004308"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_diste.rvs(mu=dm_dist,sigma=.2*dm_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import time\n",
    "num = 10000\n",
    "tmp_burm_3 = np.zeros(num)\n",
    "start = time.time()\n",
    "for i in range(int(num)):\n",
    "    tmp_burm_3[i] =tmp_diste.rvs(mu=dm_dist,sigma=.2*dm_dist)\n",
    "    print(\"\\r\",end=\"iter: {} Total time: {}\".format(i,time.time()-start),flush=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'tmp_burm' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-19-3d9e13c967a2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp_burm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp_burm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp_burm_1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp_burm_1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp_burm_2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'tmp_burm' is not defined"
     ]
    }
   ],
   "source": [
    "print(np.shape(tmp_burm))\n",
    "print(min(tmp_burm))\n",
    "print(np.shape(tmp_burm_1))\n",
    "print(min(tmp_burm_1))\n",
    "print(np.shape(tmp_burm_2))\n",
    "print(min(tmp_burm_2))\n",
    "print(np.shape(tmp_burm_3))\n",
    "print(min(tmp_burm_3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_burm_tot = np.concatenate((tmp_burm,tmp_burm_1,tmp_burm_2,tmp_burm_3))\n",
    "print(tmp_burm_tot.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#tmp_dist = NE2001DMDistPrior(px,dm_dist,.2*dm_dist)\n",
    "px = np.linspace(1e-3,10,100)\n",
    "plt.plot(px,px_prior(px,dm_dist,.2*dm_dist))\n",
    "plt.plot(px,scipy.stats.norm.pdf(px**-1,loc=dm_dist,scale=.2*dm_dist)/px**2)\n",
    "plt.plot(px,tmp_diste.pdf(px,mu=dm_dist,sigma=.2*dm_dist))\n",
    "plt.hist(tmp_burm_tot,density=True,bins=100);\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt('px_prior_1.txt', tmp_burm_tot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAHwCAYAAAD0N5r7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAABYlAAAWJQFJUiTwAAAfZUlEQVR4nO3df6zmVX0n8PdHh/hjKihulKy2TiHCEGuoYPyxbFQgy9KaNmjLf1KXlCasZFEjiUbrMmNi4qZuRfxRzbpKhc1mN2lctykqu3YWRdZ0M9S4xkEQuVQjLEWWEdiprXj2j+d7dbjcZ+a59/nee8+9z+uVPDn3+f44z5nznBne93C+32+11gIAAPThKVvdAAAA4BcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgI7u2ugGbqaruSXJikqUtbgoAADvbniQ/bq396lpPXKiAnuTEZzzjGSefeeaZJ291QwAA2LkOHTqUI0eOrOvcRQvoS2eeeebJBw8e3Op2AACwg51zzjm5/fbbl9ZzrjXoAADQEQEdAAA6IqADAEBHBHQAAOiIgA4AAB0R0AEAoCMCOgAAdERABwCAjgjoAADQEQEdAAA6IqADAEBHBHQAAOjI3AG9qp5bVZdX1eeq6rtVdaSqDlfVrVX1+1U182dU1VJVtSmv++dtKwAA9G7XCHVckuRPktyX5ECSv0ny/CRvTPKpJL9RVZe01tqM9R1Ocu0q2x+dv6kAANC3MQL6nUl+O8lftNZ+tryxqt6d5K+S/E4mYf3PZqzv4dbavhHaBQAA287cS1xaa3/ZWvvzo8P5sP3+JJ8Y3r5u3s8BAIBFMMYM+rH8w1D+dA3nPK2q3pTkV5I8luSbSb7SWnt87MYBAEBvNiygV9WuJL83vP3iGk49JckNK7bdU1WXtdZumfGzD07ZtXcN7QAAgE23kbdZ/ECSX0tyU2vtSzOe85kkF2QS0ncneWmSTybZk+QLVXXWBrSTGVXVqi8AAMazITPoVXVVknckuSPJpbOe11rbv2LTt5JcUVWPDvXtS/KGGeo5Z0q7DiY5e9b2AADAZht9Br2qrkzy4STfTnJea+2hEapdvtj0NSPUBQAA3Ro1oFfV25J8NJOZ7/OGO7mM4YGh3D1SfQAA0KXRAnpVvTPJh5J8I5Nw/sCxz1iTVw/l90asEwAAujNKQK+q92ZyUejBJBe01h48xrEnVNXeqjptxfaXVNXJqxz/okxm5ZPkxjHaCwAAvZr7ItGqenOS9yV5PMlXk1y1yp09llpr1w8/vyDJoST3ZnJ3lmWXJHlXVR1Ick+SR5KcluT1SZ6e5KYkH5y3vQAA0LMx7uLyq0P51CRvm3LMLUmuP049B5KckeRlmSxp2Z3k4SS3ZnJf9Btaa22+pgIAQN/mDuittX2Z3P5w1uOXkjxpin14CNFMDyICAICdaiMfVAQAAKyRgA4AAB0R0AEAoCMCOgAAdERABwCAjgjoAADQEQEdAAA6IqADAEBHxniSKDtM1ZOeI7Wu4z34FQBg7cygAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAju7a6AexcVTV1X2ttE1sCALB9mEEHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgI+6DvsCOdZ9yAAC2hhl0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0JG5A3pVPbeqLq+qz1XVd6vqSFUdrqpbq+r3q2pNn1FVL6yqT1fVD6vqJ1W1VFXXVtVz5m0rAAD0btcIdVyS5E+S3JfkQJK/SfL8JG9M8qkkv1FVl7TW2vEqqqrTktyW5HlJPp/kjiSvSPLWJBdV1bmttR+N0GYAAOjSGAH9ziS/neQvWms/W95YVe9O8ldJfieTsP5nM9T18UzC+VWttY8cVdcfJ3l7kvcnuWKENgMAQJfmXuLSWvvL1tqfHx3Oh+33J/nE8PZ1x6unqk5NcmGSpSQfW7H7miSPJbm0qnbP22YAAOjVRl8k+g9D+dMZjj1/KG9eJew/kuRrSZ6Z5FXjNQ8AAPoyxhKXVVXVriS/N7z94gynnDGUd07Zf1cmM+ynJ/nycT774JRde2doBwAAbJkNC+hJPpDk15Lc1Fr70gzHnzSUh6fsX97+7DnbRQeqatXtM1xLDACwo21IQK+qq5K8I5O7sFw6VrVDedwE11o7Z0q7DiY5e6T2AADA6EZfg15VVyb5cJJvJzmvtfbQjKcuz5CfNGX/iSuOAwCAHWfUgF5Vb0vy0STfyiSc37+G078zlKdP2f/ioZy2Rh0AALa90QJ6Vb0zyYeSfCOTcP7AGqs4MJQXrnz6aFU9K8m5SY4k+fqcTQUAgG6NEtCr6r2ZXBR6MMkFrbUHj3HsCVW1d3hq6M+11u5OcnOSPUmuXHHa/iS7k3y2tfbYGG0GAIAezX2RaFW9Ocn7kjye5KtJrlrlDh1LrbXrh59fkORQknszCeNHe0uS25JcV1UXDMe9Msl5mSxtec+87QUAgJ6NcReXXx3KpyZ525Rjbkly/fEqaq3dXVUvzyTwX5TkN5Pcl+S6JPvXcMEpAABsS3MH9NbaviT71nD8Un5xy8TV9n8/yWXztgsAALaj0W+zCAAArJ+ADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgI7u2ugFwtKpadXtrbZNbAgCwNcygAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0JFdW90ANl5VbXUTAACYkRl0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBAR0YJ6FX1u1X1kar6alX9uKpaVd24jnqWhnNXe90/RlsBAKBnYz1J9A+TnJXk0SQ/SLJ3jroOJ7l2le2PzlEnAABsC2MF9LdnEsy/m+S1SQ7MUdfDrbV9YzQKAAC2m1ECemvt54G8qsaoEgAAFtJYM+hjelpVvSnJryR5LMk3k3yltfb41jYLAAA2Xo8B/ZQkN6zYdk9VXdZau2WWCqrq4JRd86yNBwCADdfbbRY/k+SCTEL67iQvTfLJJHuSfKGqztq6pgEAwMbraga9tbZ/xaZvJbmiqh5N8o4k+5K8YYZ6zllt+zCzfvaczQQAgA3T2wz6NJ8YytdsaSsAAGCDbZeA/sBQ7t7SVgAAwAbbLgH91UP5vS1tBQAAbLBND+hVdUJV7a2q01Zsf0lVnbzK8S9K8tHh7Y2b0UYAANgqo1wkWlUXJ7l4eHvKUL66qq4ffn6wtXb18PMLkhxKcm8md2dZdkmSd1XVgST3JHkkyWlJXp/k6UluSvLBMdoLAAC9GusuLr+e5M0rtp06vJJJGL86x3YgyRlJXpbJkpbdSR5Ocmsm90W/obXWxmkuAAD0aZSA3lrbl8ktEGc5dilJrbL9liQzPYgIAAB2qu1ykSgAACwEAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOjPUkUdhQVU96tlWSxMNlAYCdxgw6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICO7NrqBsA8qmrV7a21TW4JAMA4zKADAEBHBHQAAOiIgA4AAB0R0AEAoCMCOgAAdERABwCAjgjoAADQEQEdAAA6IqADAEBHBHQAAOiIgA4AAB0R0AEAoCMCOgAAdERABwCAjgjoAADQEQEdAAA6IqADAEBHBHQAAOiIgA4AAB0R0AEAoCMCOgAAdGTXVjeA8VTVVjcBAIA5mUEHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgI+6Dzo407Z7wrbVNbgkAwNqYQQcAgI4I6AAA0JFRAnpV/W5VfaSqvlpVP66qVlU3rrOuF1bVp6vqh1X1k6paqqprq+o5Y7QVAAB6NtYa9D9MclaSR5P8IMne9VRSVacluS3J85J8PskdSV6R5K1JLqqqc1trPxqlxQAA0KGxlri8PcnpSU5M8i/nqOfjmYTzq1prF7fW3tVaOz/Jh5KckeT9c7cUAAA6NkpAb60daK3d1ea4RUZVnZrkwiRLST62Yvc1SR5LcmlV7V53QwEAoHM9XSR6/lDe3Fr72dE7WmuPJPlakmcmedVmNwwAADZLT/dBP2Mo75yy/65MZthPT/LlY1VUVQen7FrX2ngAANgsPc2gnzSUh6fsX97+7I1vCgAAbI2eZtCPZ/nRkMdd595aO2fVCiYz62eP2SgAABhTTzPoyzPkJ03Zf+KK4wAAYMfpKaB/ZyhPn7L/xUM5bY06AABsez0F9ANDeWFVPaFdVfWsJOcmOZLk65vdMAAA2CybHtCr6oSq2js8NfTnWmt3J7k5yZ4kV644bX+S3Uk+21p7bFMaCgAAW2CUi0Sr6uIkFw9vTxnKV1fV9cPPD7bWrh5+fkGSQ0nuzSSMH+0tSW5Lcl1VXTAc98ok52WytOU9Y7QXAAB6NdZdXH49yZtXbDt1eCWTMH51jqO1dndVvTzJ+5JclOQ3k9yX5Lok+1trD43UXgAA6NIoAb21ti/JvhmPXcovbpm42v7vJ7lsjHYBAMB209NFogAAsPAEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6MiurW4AbKaqWnV7a22TWwIAsDoz6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoyK6tbgD0oKpW3d5a2+SWAACLzgw6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHRHQAQCgI7u2ugGsXVVtdRMAANggZtABAKAjAjoAAHTEEhc4hmnLiVprm9wSAGBRjDaDXlUvrKpPV9UPq+onVbVUVddW1XPWUMdSVbUpr/vHaisAAPRqlBn0qjotyW1Jnpfk80nuSPKKJG9NclFVndta+9GM1R1Ocu0q2x8doakAANC1sZa4fDyTcH5Va+0jyxur6o+TvD3J+5NcMWNdD7fW9o3ULgAA2FbmXuJSVacmuTDJUpKPrdh9TZLHklxaVbvn/SwAANjpxphBP38ob26t/ezoHa21R6rqa5kE+Fcl+fIM9T2tqt6U5FcyCfffTPKV1trjI7QVAAC6NkZAP2Mo75yy/65MAvrpmS2gn5LkhhXb7qmqy1prt8zSoKo6OGXX3lnOBwCArTLGXVxOGsrDU/Yvb3/2DHV9JskFmYT03UlemuSTSfYk+UJVnbXuVgIAwDawGfdBX76R9HFvHN1a279i07eSXFFVjyZ5R5J9Sd4wQz3nrNqQycz62cc7HwAAtsoYM+jLM+QnTdl/4orj1uMTQ/maOeoAAIDujRHQvzOUp0/Z/+KhnLZGfRYPDKU7wQAAsKONEdAPDOWFVfWE+qrqWUnOTXIkydfn+IxXD+X35qgDAAC6N3dAb63dneTmTC7kvHLF7v2ZzHp/trX2WJJU1QlVtXd4+ujPVdVLqurklfVX1YuSfHR4e+O87QUAgJ6NdZHoW5LcluS6qrogyaEkr0xyXiZLW95z1LEvGPbfm0moX3ZJkndV1YEk9yR5JMlpSV6f5OlJbkrywZHaCwAAXRoloLfW7q6qlyd5X5KLkvxmkvuSXJdkf2vtoRmqOZDJPdVflsmSlt1JHk5yayb3Rb+htXbcO8EAAMB2NtptFltr309y2QzHLeUXt148evstSWZ6EBEAAOxUY1wkCgAAjERABwCAjgjoAADQEQEdAAA6IqADAEBHRruLCyySqifdiChJ4k6gAMC8zKADAEBHBHQAAOiIgA4AAB0R0AEAoCMCOgAAdERABwCAjgjoAADQEQEdAAA6IqADAEBHBHQAAOiIgA4AAB0R0AEAoCO7troBsJNU1arbW2ub3BIAYLsygw4AAB0R0AEAoCMCOgAAdERABwCAjgjoAADQEQEdAAA6IqADAEBHBHQAAOiIgA4AAB0R0AEAoCMCOgAAdGTXVjcAFkFVTd3XWtvElgAAvTODDgAAHRHQAQCgIwI6AAB0xBr0jh1r3TI7x7Tv2dp0AFhMZtABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiPugQ6fcHx0AFpMZdAAA6IiADgAAHRHQAQCgIwI6AAB0REAHAICOCOgAANARAR0AADoioAMAQEcEdAAA6IiADgAAHdm11Q0A1qaqVt3eWtvklgAAG8EMOgAAdMQMOuwQZtYBYGcwgw4AAB0R0AEAoCOWuMAOZ+kLAGwvZtABAKAjAjoAAHTEEhdYUJa+AECfzKADAEBHzKADT2BmHQC2loDegWmBCACAxWOJCwAAdMQMOjCTtf6fHktiAGB9zKADAEBHRgvoVfXCqvp0Vf2wqn5SVUtVdW1VPWcr6gG2VlWN9gKARTLKEpeqOi3JbUmel+TzSe5I8ookb01yUVWd21r70WbVAwAA29VYM+gfzyRUX9Vau7i19q7W2vlJPpTkjCTv3+R6gAVgxh2AnWjugF5Vpya5MMlSko+t2H1NkseSXFpVuzejHmDnEcQBWCRjzKCfP5Q3t9Z+dvSO1tojSb6W5JlJXrVJ9XRLyIDNMeb6d2voAdhsY6xBP2Mo75yy/65MZsZPT/LlTagnVXVwyq6zDh06lHPOOedYpwOMSkjv29lnn72m42+//fYNrX+tjtWeaZ897ZyNbut2sl36aLu0czP01heHDh1Kkj3rOXeMgH7SUB6esn95+7M3qZ5jefzIkSOHb7/99qU56hjD3qG8Y0tbsb3pw3Hox3Hox/ltWR+uNXD3Vv8KT+jHtX72Jre1V8cci9uljzpoZzf/Lm5hX+xJ8uP1nLgZDypanjqa96klM9fTWut6inx5hr/3dvZMH45DP45DP85PH45DP85PH45DP85njDXoyzPbJ03Zf+KK4za6HgAA2LbGCOjfGcrTp+x/8VBOW1s+dj0AALBtjRHQDwzlhVX1hPqq6llJzk1yJMnXN6keAADYtuYO6K21u5PcnMlC+CtX7N6fZHeSz7bWHkuSqjqhqvYOTw1ddz0AALATjXWR6FuS3Jbkuqq6IMmhJK9Mcl4mS1Lec9SxLxj235sn33pmLfUAAMCOU63Ne3OVoaKqX07yviQXJXlukvuS/Jck+1trDx113J4k9yS5t7W2Z731AADATjRaQAcAAOY3xkWiAADASAR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6COoqhdW1aer6odV9ZOqWqqqa6vqOWuoY6mq2pTX/RvZ/h5U1e9W1Ueq6qtV9ePhz33jOuua+/vYjsbqw0Uei1X13Kq6vKo+V1XfraojVXW4qm6tqt+vqjX9m7nAY3G0flzw8fhvqurLVfX9oQ8fqqq/rqprquq5a6xrIcdiMl4/LvJYXE1VXXrUn//yNZ67sONxVu6DPqeqOi2Tp58+L8nnk9yR5BWZPP30O0nOba39aIZ6lpI8O8m1q+x+tLX2wXFa3Keq+kaSs5I8muQHSfYm+Q+ttTetsZ5Rvo/taMQ+XMqCjsWquiLJn2TygLQDSf4myfOTvDHJSUn+LMklbYZ/OBd8LI7Zj0tZ3PH490luT/LtJA8k2Z3kVUlenuSHSV7VWvv+DPUs7FhMRu3HpSzoWFypJg+V/N9Jnprkl5L8QWvtUzOeu9DjcWatNa85Xkm+lKQl+Vcrtv/xsP0TM9azlGRpq/88W9iP5yV5cZJK8rqh727cqu9jO75G7MOFHYtJzk/yW0mesmL7KZmEzJbkd2asa5HH4pj9uMjj8elTtr9/6MOPz1jPwo7FkftxYcfiin6oJP89yd1J/mjow8vXcP5Cj8dZX5a4zKGqTk1yYSZ/aT+2Yvc1SR5LcmlV7d7kpm07rbUDrbW72vC3dD0W/fsYow8XXWvtL1trf95a+9mK7fcn+cTw9nXHq8dYHKcfF11r7e+m7PrPQ/ni49Wx6GMxGacfeYKrMvkl/LJMxs/MjMfZCejzOX8ob17lP0SPJPlakmdm8r/SZvG0qnpTVb27qt5aVedV1VNHbO9ON/b3sciMxSf7h6H86QzHGovTraUflxmPT/RbQ/nNGY41FqdbSz8uW+ixWFVnJvlAkg+31r6yjiqMxxnt2uoGbHNnDOWdU/bflclviqcn+fIM9Z2S5IYV2+6pqstaa7esr4kLZezvY5EZi0epql1Jfm94+8UZTjEWV7GOfly20OOxqq7OZJ3vSZmsm/6nmYTKD8xwurE4mLMfly3sWBz+/t6QyTK1d6+zGuNxRmbQ53PSUB6esn95+7NnqOszSS7I5C//7iQvTfLJJHuSfKGqzlp3KxfHmN/HIjMWn+wDSX4tyU2ttS/NcLyxuLq19mNiPCbJ1Zn87/+3ZRIqv5jkwtba385wrrH4C/P0Y2Is/uskL0vyL1prR9ZZh/E4IwF9Y9VQHndNcGtt/7Bu8/+01v5fa+1brbUrMrlo4hlJ9m1gOxfFzN/HIjMWn6iqrkryjkzuNHDpWNUO5cKMxfX2o/GYtNZOaa1VJsHwjUlOTfLXVXX2CNUvzFictx8XeSxW1SsymTX/t621/7mRHzWUO348Ho+APp/l3/ROmrL/xBXHrcfyBVWvmaOORbEZ38ciW7ixWFVXJvlwJrdnO6+19tCMpxqLR5mjH49l4cbjEAw/l8kSgOcm+ewMpxmLK6yzH49lR4/Fo5a23JnkvXNWZzzOSECfz3eG8vQp+5evDJ+21moWDwzlwl/RPIPN+D4W2UKNxap6W5KPJvlWJqFyLQ8iMRYHc/bjsSzUeDxaa+3eTH7ZeUlV/aPjHG4sTrHGfjyWnT4WfymT8XNmkr87+gFNmSwZSpJ/N2y79jh1GY8zcpHofA4M5YVV9ZSjr0iuqmclOTfJkSRfn+MzXj2U35ujjkWxGd/HIluYsVhV78xkvfQ3kvyz1tqDa6zCWMwo/XgsCzMep/jHQ/n4cY4zFo9t1n48lp0+Fn+S5N9P2Xd2JuvSb80kfB9v+YvxOCMz6HNord2d5OZMLhC5csXu/Zn8Nv3Z1tpjSVJVJ1TV3uEpWj9XVS+pqpNX1l9VL8pk5ilJ1vXY+51oWj+u9ftYZMbidFX13kxC5cEkFxwrVBqL043Rj4s8Hof+OGWV7U+pqvdn8hTG21pr/3fYbiyuYqx+XOSx2Fo70lq7fLVXkv86HPanw7b/lBiPYyjPNJnPKo+sPZTklZk81fHOJP+kDY+srao9Se5Jcm9rbc9RdexL8q5MfrO8J8kjSU5L8vokT09yU5I3tNb+fjP+TFuhqi5OcvHw9pQk/zyT2YivDtsebK1dPRy7J6v047Bv5u9jpxmjDxd9LFbVm5Ncn8ls2key+jrIpdba9cPxe2IsPslY/bjI43FYGvRHSb6SyRMbf5Tk+Ulem8nFjfdn8ovPt4fj98RYfJKx+nGRx+KxDP1yTZI/aK196qjte2I8zqd18DjT7f5K8suZ3H7pviR/n+TeTC6IOnnFcXsyuTJ5acX21yb5j5nc3eDhTB7k8bdJ/lsm9wyurf4zbkIf7hv6Ztpr6ahjV+3HtX4fO+01Rh8u+licoQ9bkv9hLG5OPy7yeMzkdpQfy2R50IOZPNjpcJL/NfTvTP99MRbH6cdFHovH6d/lv+uXz9KPR+1fyPG4lpcZdAAA6Ig16AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0BEBHQAAOiKgAwBARwR0AADoiIAOAAAdEdABAKAjAjoAAHREQAcAgI4I6AAA0JH/D/BgUftr9KwdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 248,
       "width": 372
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = np.loadtxt('px_prior_1.txt')\n",
    "plt.hist(a,density=True,bins=100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NE2001DMDist_rv_gen(rv_continuous):\n",
    "    \"Gaussian distribution\"\n",
    "    def _pdf(self,x, mu, sigma):\n",
    "        return (1/x**2)*scipy.stats.norm.pdf(x**-1,loc=mu,scale=sigma)\n",
    "    \n",
    "def NE2001DMDistPrior(value, mu=1., sigma=.2):\n",
    "    \"\"\"Prior function for DE2001 DMDistance parameters.\"\"\"\n",
    "    rv_dist = NE2001DMDist_gen(shapes='mu,sigma')\n",
    "    return rv_dist.pdf(value,mu,sigma)\n",
    "\n",
    "\n",
    "def NE2001DMDistSampler(mu=1., sigma=.2, size=None):\n",
    "    \"\"\"Sampling function for DE2001 DMDistance parameters.\"\"\"\n",
    "    rv_dist = NE2001DMDist_gen(shapes='mu,sigma')\n",
    "    return rv_dist.rvs(mu=mu,sigma=sigma,size=size)\n",
    "\n",
    "def NE2001DMDist(mu=1., sigma=.2, size=None):\n",
    "    \"\"\"Class factory for bounded Normal parameters.\"\"\"\n",
    "\n",
    "    class NE2001DMDist(parameter.Parameter):\n",
    "        _prior = parameter.Function(\n",
    "            NE2001DMDistPrior, mu=mu, sigma=sigma,\n",
    "        )\n",
    "        _sampler = staticmethod(NE2001DMDistSampler)\n",
    "        _size = size\n",
    "        _mu = mu\n",
    "        _sigma = sigma\n",
    "        _pmin = pmin\n",
    "        _pmax = pmax\n",
    "\n",
    "        def __repr__(self):\n",
    "            return \"{}: NE2001DMDist({},{}, [{},{}])\".format(\n",
    "                self.name, mu, sigma, pmin, pmax\n",
    "            ) + (\"\" if self._size is None else \"[{}]\".format(self._size))\n",
    "\n",
    "    return NE2001DMDist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}