distribution

gistfile1.py

#!/usr/bin/env python
from matplotlib import pyplot as plt
from exparser import TraceKit as tk
import numpy as np

def findPeak(a, searchRange, smoothing=201, N=1000, plot=False):

    """
    desc:
        Get the peak of a series of data.

    arguments:
        a:
            desc:   The data as a 1D array.
            type:   ndarray
        searchRange:
            desc:   The range to search in as a (min, max) tuple.
            type:   tuple
        smoothing:
            desc:   A smoothing value in case the data is rough.
            type:   int
        N:
            desc:   The resolution of the density plot used to find the peak.
            type:   N
        plot:
            desc:   Indicates whether a debug plot should be shown.
            type:   bool

    returns:
        desc:   The peak of the distribution.
        type:   float
    """

    # Create a cumulative distribution
    a.sort()
    ax = np.linspace(searchRange[0], searchRange[1], N)
    ay = np.empty(N)
    for i, x in enumerate(ax):
        y = np.sum(a <= x)
        ay[i] = y
    # Smooth the cumulative distribution
    asy = tk.smooth(ay, windowLen=smoothing)
    # Get the derivative of the cumulative distribution
    adx = .5*(ax[1:]+ax[:-1])
    ady = np.abs(asy[1:]-asy[:-1])
    # Get the peak (or through)
    iMax = np.where(ady == np.max(ady))[0]
    xMax = adx[iMax]
    if plot:
        plt.subplot(2,1,1)
        plt.plot(ax, ay, color='green')
        # Plot the smoothed
        plt.plot(ax, asy, color='blue')
        plt.subplot(2,1,2)
        plt.plot(adx, ady, color='red')
        plt.axvline(xMax)
        plt.show()
    return xMax

# Generate data
a = np.random.normal(0, .5, size=100)
peak = findPeak(a, (-2, 2), plot=True)