Skip to content

Commit

Permalink
Merge pull request #301 from sblunt/add_mapfit
Browse files Browse the repository at this point in the history 
Adding scipy.optimize.minimize to modify initial positions
  • Loading branch information
@sblunt
sblunt committed Nov 2, 2021
2 parents 665ca48 + 7cb0f09 commit b5b2fcc
Showing 1 changed file with 119 additions and 68 deletions.
187 changes: 119 additions & 68 deletions docs/tutorials/Modifying_MCMC_initial_positions.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
"source": [
"# Modifying MCMC Initial Positions\n",
"\n",
"by Henry Ngo (2019) & Sarah Blunt (2021)\n",
"by Henry Ngo (2019) & Sarah Blunt (2021) & Mireya Arora (2021)\n",
"\n",
"When you set up the MCMC Sampler, the initial position of your walkers are randomly determined. Specifically, they are uniformly distributed in your Prior phase space. This tutorial will show you how to change this default behaviour so that the walkers can begin at locations you specify. For instance, if you have an initial guess for the best fitting orbit and want to use MCMC to explore posterior space around this peak, you may want to start your walkers at positions centered around this peak and distributed according to an N-dimensional Gaussian distribution. \n",
"\n",
Expand All @@ -30,12 +30,12 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"from scipy.optimize import minimize as mn\n",
"import orbitize\n",
"from orbitize import driver\n",
"import multiprocessing as mp"
Expand All @@ -52,7 +52,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -92,7 +92,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 32,
"metadata": {},
"outputs": [
{
Expand All @@ -118,14 +118,14 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['__abstractmethods__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__slots__', '__str__', '__subclasshook__', '__weakref__', '_abc_impl', '_fill_in_fixed_params', '_logl', '_update_chains_from_sampler', 'check_prior_support', 'chop_chains', 'curr_pos', 'custom_lnlike', 'examine_chains', 'fixed_params', 'has_corr', 'lnlike', 'num_params', 'num_temps', 'num_threads', 'num_walkers', 'priors', 'results', 'run_sampler', 'sampled_param_idx', 'system', 'use_c', 'use_gpu', 'use_pt', 'validate_xyz_positions']\n"
"['__abstractmethods__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__slotnames__', '__slots__', '__str__', '__subclasshook__', '__weakref__', '_abc_impl', '_fill_in_fixed_params', '_logl', '_update_chains_from_sampler', 'check_prior_support', 'chop_chains', 'curr_pos', 'custom_lnlike', 'examine_chains', 'fixed_params', 'has_corr', 'lnlike', 'num_params', 'num_temps', 'num_threads', 'num_walkers', 'priors', 'results', 'run_sampler', 'sampled_param_idx', 'system', 'use_pt', 'validate_xyz_positions']\n"
]
}
],
Expand All @@ -151,7 +151,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 34,
"metadata": {},
"outputs": [
{
Expand All @@ -160,7 +160,7 @@
"(5, 30, 8)"
]
},
"execution_count": 5,
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -178,23 +178,23 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[2.81112329e+01 6.78182187e-01 2.01780023e+00 5.72736907e+00\n",
" 4.17834790e+00 4.47381587e-01 5.15315635e+01 1.74688938e+00]\n",
" [4.15746056e-01 4.01910619e-01 1.86025997e+00 2.87105502e-01\n",
" 5.95428436e+00 6.12526807e-01 5.13313374e+01 1.82416751e+00]\n",
" [1.22068925e+02 4.67547720e-03 1.68530675e+00 8.42270456e-01\n",
" 3.57627204e+00 5.84369274e-01 5.13120124e+01 1.77019622e+00]\n",
" [2.30625205e+02 7.42423534e-01 1.26033986e+00 8.05090822e-01\n",
" 5.56480261e+00 4.29398719e-01 5.15298178e+01 1.75015701e+00]\n",
" [3.58583894e-03 1.04218960e-01 4.05009363e-01 5.72259383e+00\n",
" 3.12529358e+00 5.89550947e-01 5.16473143e+01 1.68496551e+00]]\n"
"[[2.96649493e-03 6.46389191e-02 2.05568393e+00 5.59133817e+00\n",
" 6.15061921e+00 9.04132185e-01 5.13464691e+01 1.71444962e+00]\n",
" [8.45997206e+03 2.28786525e-01 1.05496996e+00 5.40785279e-01\n",
" 3.12679020e+00 9.83935240e-01 5.15802933e+01 1.79068792e+00]\n",
" [2.06674323e-01 4.14942416e-01 9.10528634e-01 5.12112580e+00\n",
" 1.40423713e+00 2.11681325e-01 5.13803625e+01 1.67395370e+00]\n",
" [7.58805546e-02 2.95768231e-01 1.24678730e+00 1.77255326e+00\n",
" 3.11043407e+00 6.94950535e-01 5.14178357e+01 1.83083582e+00]\n",
" [6.64985800e-02 5.30360650e-01 2.72837627e+00 3.58772106e+00\n",
" 4.25024656e+00 2.51044429e-01 5.14502371e+01 1.80271446e+00]]\n"
]
}
],
Expand Down Expand Up @@ -235,7 +235,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 36,
"metadata": {},
"outputs": [
{
Expand All @@ -252,7 +252,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -285,7 +285,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -299,106 +299,157 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.2) Validate your new positions\n",
"\n",
"Drawing from a normal distribution can cause your walkers to start outside of your prior space. See the Modifying Priors tutorial for information on how to interact with the prior objects, which would allow you to find the limits on each parameter set by the priors etc. \n",
"\n",
"Here, let's do something more simple and just check that all values are physically valid. In this tutorial, eccentricity is the most likely problem because the distribution from Blunt et al. 2017 was very non-Gaussian but we are applying a Gaussian distribution with centre at 0.0151 but a spread of 0.175, so almost half of the generated numbers will be negative! \n",
"\n",
"So, let's keep the default eccentricity values of the walkers (originally uniformly distributed from 0 to 1). You can make a copy of the current position into the `new_pos` array. Another option (not shown here) would be to generate a different distribution (e.g. Poisson) for this parameter instead."
"### 3.2) Using an optimizer to obtain a best fit value"
]
},
{
"cell_type": "code",
"execution_count": 10,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"ecc_ind = my_driver.system.param_idx['ecc1']\n",
"new_pos[:,:,ecc_ind] = np.copy(my_driver.sampler.curr_pos[:,:,ecc_ind])"
"Other optimizing software can also be used to generate intial positions. Depending on the quality of data collected and whether a suitable guess array of parameters can be made, different optimizing software can provide better best fit values for for MCMC walkers. Below you will find a few options that cater to different scenarios."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Randomizing some angles\n",
"You could also just change some values in the `new_pos` arrays. For instance, the `aop` and `pan` angles are degenerate by 180 degrees (i.e. 30 degrees is the same as 210 degrees). If you are getting values from a previous fit, they might have already been clustered around one of the values. So, we can randomly take about half of the values and add or subtract 180 degrees."
"### 3.2a) Using `scipy.optimize.minimize`"
]
},
{
"cell_type": "code",
"execution_count": 11,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"# Find indices of the two angles\n",
"aop_ind = my_driver.system.param_idx['aop1']\n",
"pan_ind = my_driver.system.param_idx['pan1']\n",
"Assuming the data obtained allows for a suitable guess to be made for each parameter, a scipy.optimize.minimize software can be used to generate a best fit value. You may want to skip this step and input your guess values directly into MCMC's initial walker positions, however scipy can help refine the guess parameters.\n",
"\n",
"# Get the shape of curr_pos without the last dimension (the list of parameters)\n",
"select_arr_shape = curr_pos_shape[:-1] # (n_temps,n_walkers)\n",
"First, we define a new log liklihood function function `neg_logl` based on the guess values we have. \n",
"Note, since we have predefined a good guess, from the aforementioned Table, as `walker_centres` we will continue to use it as a guess array for examples below."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[4.90426906e+01 9.99455386e-06 2.48501010e+00 1.60085022e+00\n",
" 2.33034111e+00 7.69934177e-01 5.14404237e+01 1.74922411e+00]\n"
]
}
],
"source": [
"# The following code performs a minimization whereas the log liklihood function is based on maximization so we redefine the\n",
"# likelihood function is redefined to return -x to make this a minization scenario \n",
"\n",
"# Draw a random number from 0 to 1, and mark index for replacement if > 0.5\n",
"replace_index = np.random.uniform(size=select_arr_shape)>0.5\n",
"m = my_driver.sampler\n",
"\n",
"# Replace aop values selected with current values plus 180 degrees\n",
"new_pos[replace_index,aop_ind] = new_pos[replace_index,aop_ind] + 180.0\n",
"def neg_logl(paramarray):\n",
" x = m._logl(paramarray, include_logp=True) #set include_logp to true to include guess array in likelihood function\n",
" \n",
" return -x \n",
"\n",
"# This may cause values to be larger than 360; find these and subtract 360\n",
"wrap_ind = new_pos[:,:,aop_ind] > 360\n",
"new_pos[wrap_ind] = new_pos[wrap_ind] - 360\n",
"\n",
"# Repeat all of the above for the pan angle\n",
"replace_index = np.random.uniform(size=select_arr_shape)>0.5\n",
"new_pos[replace_index,pan_ind] = new_pos[replace_index,pan_ind] + 180.0\n",
"wrap_ind = new_pos[:,:,pan_ind] > 360\n",
"new_pos[wrap_ind] = new_pos[wrap_ind] - 360"
"guessarray = walker_centres\n",
"results = mn(neg_logl, guessarray, method=\"Powell\")\n",
"print(results.x) #results.x is the best fit value"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Additional checks\n",
"\n",
"The following function can be used to identify walkers that have been initialized outside of the appropriate prior probability space. It will raise a `ValueError` if walkers are initialized outside of the priors. You should update your positions until this method runs without raising an error."
"In our trials, Powell has given the best results, but you may replace it with a different minimizing method depending on your need. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"cell_type": "markdown",
"metadata": {},
"outputs": [],
"source": [
"my_driver.sampler.check_prior_support()"
"### 3.3) Scattering walkers"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We should continue investigating which parameters are being initialized outside of the prior space until this function returns empty lists."
"To set up MCMC so that it explores the nearby probablity space thoroughly and finds the global minimum, you can scatter the initial positions of the walkers around the best fit value. This can be done by adding random numbers to `results.x`\n",
"\n",
"This section overrides walker_1sigmas and creates a spread of `new_pos` in a different manner than above.\n",
"The following is a template based on the aforementioned Table. The scatter is created using a variety of methods, we recommend reviewing the code to ensure it is compatible to your data."
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"new_pos = np.random.standard_normal(curr_pos_shape)*0.03 + results.x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.3) Update sampler position\n",
"### 3.4) Update sampler position\n",
"\n",
"After generating and validating your new walker positions, through whatever methods you choose, it's now time to update the sampler object to have its `curr_pos` be your new positions."
]
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"my_driver.sampler.curr_pos = np.copy(new_pos)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.5) Validate your new positions\n",
"\n",
"Drawing from a normal distribution can cause your walkers to start outside of your prior space. See the Modifying Priors tutorial for information on how to interact with the prior objects, which would allow you to find the limits on each parameter set by the priors etc. \n",
"\n",
"Here, let's do something more simple and just check that all values are physically valid. After this we can begin to correct them.\n",
"\n",
"The following function can be used to identify walkers that have been initialized outside of the appropriate prior probability space. It will raise a `ValueError` if walkers are initialized outside of the priors. You should update your positions until this method runs without raising an error."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "Attempting to start with walkers outside of prior support: check parameter(s) 1",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-42-1b66a275c5b0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmy_driver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msampler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcheck_prior_support\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[0;32m~/Projects/orbitize/orbitize/sampler.py\u001b[0m in \u001b[0;36mcheck_prior_support\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 1094\u001b[0m \u001b[0;31m# Throw our ValueError if necessary,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1095\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbad_parameters\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1096\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Attempting to start with walkers outside of prior support: check parameter(s) \"\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m', '\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbad_parameters\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 1097\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1098\u001b[0m \u001b[0;31m# We're not done yet, however. There may be errors in covariant priors; run a check for that.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: Attempting to start with walkers outside of prior support: check parameter(s) 1"
]
}
],
"source": [
"my_driver.sampler.check_prior_support()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We should continue investigating which parameters are being initialized outside of the prior space until this function returns empty lists."
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -412,7 +463,7 @@
"hash": "e899b22145868d3cd465733d82c36c2ae3ac0d3591d6a0807ec2e5e577a9cf5c"
},
"kernelspec": {
"display_name": "Python 3.7.5 64-bit ('python3.7': conda)",
"display_name": "Python 3.7.5 64-bit",
"name": "python3"
},
"language_info": {
Expand Down

0 comments on commit b5b2fcc

Please sign in to comment.