Priors for Bayesian analysis

Astromodels supports the definition of priors for all parameters in your model. You can use as prior any function (although of course not all functions should be used this way, but the choice is up to you).

First let's define a simple model containing one point source (see the "Model tutorial" for more info):



In [1]:

    
from astromodels import *

# Create a point source named "pts1"
pts1 = PointSource('pts1',ra=125.23, dec=17.98, spectral_shape=powerlaw())

# Create the model
my_model = Model(pts1)









    



/home/giacomov/software/canopy-env/lib/python2.7/site-packages/astromodels-0.1-py2.7.egg/astromodels/functions/functions.py:39: NaimaNotAvailable: The naima package is not available. Models that depend on it will not be available

Now let's assign uniform priors to the parameters of the powerlaw function. The function uniform_prior is defined like this:



In [2]:

    
uniform_prior.info()









    






description: A function which is constant on the interval lower_bound - upper_bound and 0 outside the interval. The extremes of the interval are counted as part of the interval.

formula: $ f(x)=\begin{cases}0 & x < \text{lower_bound} \\\text{value} & \text{lower_bound} \le x \le \text{upper_bound} \\ 0 & x > \text{upper_bound} \end{cases}$

default parameters: 


lower_bound: 


value: 0.0

desc: Lower bound for the interval

min_value: None

max_value: None

unit: 

delta: 0.1

free: True





upper_bound: 


value: 1.0

desc: Upper bound for the interval

min_value: None

max_value: None

unit: 

delta: 0.3

free: True





value: 


value: 1.0

desc: Value in the interval

min_value: None

max_value: None

unit: 

delta: 0.3

free: True

We can use it as such:



In [12]:

    
# Set 'lower_bound' to -10, 'upper bound' to 10, and leave the 'value' parameter 
# to the default value
pts1.spectrum.main.powerlaw.K.prior = log_uniform_prior(lower_bound = 1e-15, upper_bound=1e-7)

# Display it
pts1.spectrum.main.powerlaw.K.display()

# Set 'lower_bound' to -10, 'upper bound' to 0, and leave the 'value' parameter 
# to the default value
pts1.spectrum.main.powerlaw.index.prior = uniform_prior(lower_bound = -10, upper_bound=0)

pts1.spectrum.main.powerlaw.index.display()









    




Parameter K = 1.0 [1 / (cm2 keV s)]
(min_value = None, max_value = None, delta = 0.3, free = True) [prior: log_uniform_prior]






    




Parameter index = -2.0 []
(min_value = -10.0, max_value = 10.0, delta = 0.6, free = True) [prior: uniform_prior]

Now we can evaluate the prior simply as:



In [13]:

    
# Create a short cut to avoid writing too much
po = pts1.spectrum.main.powerlaw

# Evaluate the prior in 2.3e-5
point = 2.3e-21
prior_value1 = po.K.prior(point * po.K.unit)

# Equivalently we can use the fast call with no units
prior_value2 = po.K.prior.fast_call(point)

assert prior_value1 == prior_value2

print("The prior for logK evaluate to %s in %s" % (prior_value1, point))









    



The prior for logK evaluate to 0.0 in 2.3e-21

Let's plot the value of the prior at some random locations:



In [14]:

    
# You need matplotlib installed for this
import matplotlib.pyplot as plt

# This is for the IPython notebook
%matplotlib inline

# Let's get 500 points uniformly distributed between -20 and 20

random_points = np.logspace(-30,2,50)

plt.loglog(random_points,pts1.spectrum.main.powerlaw.K.prior.fast_call(random_points), '.' )
#plt.xscale("log")
#plt.ylim([-0.1,1.2])
plt.xlabel("value of K")
plt.ylabel("Prior")









    Out[14]:





<matplotlib.text.Text at 0x6566ad0>



In [ ]: