def gaussian_membership(bin_center, bin_width, smoothness):
    r"""Create a gaussian membership function for a fuzzy histogram bin.
    
    Parameters
    ----------
    bin_center : number
        The center of the bin of which to compute the membership function.
    bin_width : number
        The width of a single bin (all expected to be of equal width).
    smoothness : number, optional
        The smoothness of the function; determines the neighbourhood affected.
        See below and `fuzzy_histogram` for a more detailed explanation
        
    Returns
    -------
    gaussian_membership : function
        The cumulative density function of the desired gaussian.

    Notes
    -----
    Since the gaussian membership function is infinite, it is
    not actually true that it does not contribute to bins outside of the neighbourhood
    range. But the contribution is so marginal (:math:`eps <= 0.001` per value) that it can be
    safely ignored.
    
    The gaussian membership function is defined as
    
    .. math::
    
        \mu_{gauss}(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{(x-\zeta)^2}{2\sigma^2}}

    Since the gaussian distributions can not be formed to sum up to one at each point of
    the x-axis, their cumulative density functions (CDF) are used instead. For more
    details on CDF see http://en.wikipedia.org/wiki/Normal_distribution .
    
    The gaussian and therefore the CDF are centered above the requested value instead of
    the bin center. Then the CDF value for the left side of the bin is subtracted from
    the CDF value returned for the right side. The result is the integral under the
    gaussian with :math:`\mu/\zeta = value` with the bin-sides as the integral borders.
    
    This approach might seem a little bit unintuitive, but is the best possible for
    gaussian membership functions. The following graph gives a graphical example of the
    computation of each values bin membership
    
    .. .. image:: images/gaussian_01.png 
    
    "Trapezium functions (1)"
    
    where the bin_width is 1, one bin between each of the x tics (e.g. [-1, 0], [0, 1],
    etc.). The value which membership should be computed is marked by a yellow bar at
    :math:`x = 0.3`. Its membership in each bin is defined by the integral under the gaussian
    centered at the value (i.e. 0.3). The purple area therefore defines its membership in
    the [-2,-1] bin, the red area its membership in the [-1,0] bin, etc.
    Since the gaussian is guaranteed to have an infinite integral of 1, the some of the
    contributions of a value to all bins is one.

    For computation the function normalizes all values to a bin_width of 1, which can
    introduce marginal rounding errors.
    """
    if smoothness > 10 or smoothness < 1./10: raise AttributeError('the gaussian membership function supports only smoothnesses between 1/10 and 5.')
    
    bin_width = float(bin_width)
    bin_center = bin_center / bin_width
    start = bin_center - 0.5
    end = bin_center + 0.5
    sigma = _gaussian_membership_sigma(smoothness)

    def fun(x):
        return scipy.stats.norm.cdf(end, x / bin_width, sigma) - scipy.stats.norm.cdf(start, x / bin_width, sigma) # x, mu, sigma
    
    return fun