This documentation is for astroML version 0.2

This page


astroML Mailing List

GitHub Issue Tracker


Scipy 2012 (15 minute talk)

Scipy 2013 (20 minute talk)


If you use the software, please consider citing astroML. astroML.density_estimation.bayesian_blocks

astroML.density_estimation.bayesian_blocks(t, x=None, sigma=None, fitness='events', **kwargs)

Bayesian Blocks Implementation

This is a flexible implementation of the Bayesian Blocks algorithm described in Scargle 2012 [R14]

Parameters :

t : array_like

data times (one dimensional, length N)

x : array_like (optional)

data values

sigma : array_like or float (optional)

data errors

fitness : str or object

the fitness function to use. If a string, the following options are supported:

  • ‘events’ : binned or unbinned event data

    extra arguments are p0, which gives the false alarm probability to compute the prior, or gamma which gives the slope of the prior on the number of bins.

  • ‘regular_events’ : non-overlapping events measured at multiples

    of a fundamental tick rate, dt, which must be specified as an additional argument. The prior can be specified through gamma, which gives the slope of the prior on the number of bins.

  • ‘measures’ : fitness for a measured sequence with Gaussian errors

    The prior can be specified using gamma, which gives the slope of the prior on the number of bins. If gamma is not specified, then a simulation-derived prior will be used.

Alternatively, the fitness can be a user-specified object of type derived from the FitnessFunc class.

Returns :

edges : ndarray

array containing the (N+1) bin edges

See also

histogram plotting function which can make use of bayesian blocks.


[R14](1, 2) Scargle, J et al. (2012)


Event data:

>>> t = np.random.normal(size=100)
>>> bins = bayesian_blocks(t, fitness='events', p0=0.01)

Event data with repeats:

>>> t = np.random.normal(size=100)
>>> t[80:] = t[:20]
>>> bins = bayesian_blocks(t, fitness='events', p0=0.01)

Regular event data:

>>> dt = 0.01
>>> t = dt * np.arange(1000)
>>> x = np.zeros(len(t))
>>> x[np.random.randint(0, len(t), len(t) / 10)] = 1
>>> bins = bayesian_blocks(t, fitness='regular_events', dt=dt, gamma=0.9)

Measured point data with errors:

>>> t = 100 * np.random.random(100)
>>> x = np.exp(-0.5 * (t - 50) ** 2)
>>> sigma = 0.1
>>> x_obs = np.random.normal(x, sigma)
>>> bins = bayesian_blocks(t, fitness='measures')