Creating gaussians of fixed width and std

Question

I am trying to make every point above 25.2 a Gaussian peak with the width of 2 on the x axis. enter image description here

not so sure how to generate the Gaussian curves in python.

Answer 1

Full example of how to generate a Gaussian distribution, for an arbitrary number of axis and number of center locations. It requires the packages matplotlib, scipy and numpy.

The module can be controlled by:

• dim for the number of dimensions (axis).
• fwhm full width half maximum (estimates the width of the Gaussian distribution.)
• centers a np.array or list of the indices, that are the center(s) of the Gaussian distribution.

import matplotlib.cm as mpl_cm
import matplotlib.colors as mpl_colors
import matplotlib.pyplot as plt
import numpy as np

from scipy.spatial.distance import cdist


class Gaussian:
    def __init__(self, size):
        self.size = size
        self.center = np.array(self.size) / 2
        self.axis = self._calculate_axis()

    def _calculate_axis(self):
        """
            Generate a list of rows, columns over multiple axis.

            Example:
                Input: size=(5, 3)
                Output: [array([0, 1, 2, 3, 4]), array([[0], [1], [2]])]
        """
        axis = [np.arange(size).reshape(-1, *np.ones(idx, dtype=np.uint8))
                for idx, size in enumerate(self.size)]
        return axis

    def update_size(self, size):
        """ Update the size and calculate new centers and axis.  """
        self.size = size
        self.center = np.array(self.size) / 2
        self.axis = self._calculate_axis()

    def create(self, dim=1, fwhm=3, center=None):
        """ Generate a gaussian distribution on the center of a certain width.  """
        center = center if center is not None else self.center[:dim]
        distance = sum((ax - ax_center) ** 2 for ax_center, ax in zip(center, self.axis))
        distribution = np.exp(-4 * np.log(2) * distance / fwhm ** 2)
        return distribution

    def creates(self, dim=2, fwhm=3, centers: np.ndarray = (None,)):
        """ Combines multiple gaussian distributions based on multiple centers.  """
        centers = np.array(centers).T
        indices = np.indices(self.size).reshape(dim, -1).T

        distance = np.min(cdist(indices, centers, metric='euclidean'), axis=1)
        distance = np.power(distance.reshape(self.size), 2)

        distribution = np.exp(-4 * np.log(2) * distance / fwhm ** 2)
        return distribution

    @staticmethod
    def plot(distribution, show=True):
        """ Plotter, in case you do not know the dimensions of your distribution, or want the same interface.  """
        if len(distribution.shape) == 1:
            return Gaussian.plot1d(distribution, show)
        if len(distribution.shape) == 2:
            return Gaussian.plot2d(distribution, show)
        if len(distribution.shape) == 3:
            return Gaussian.plot3d(distribution, show)
        raise ValueError(f"Trying to plot {len(distribution.shape)}-dimensional data, "
                         f"Only 1D, 2D, and 3D distributions are valid.")

    @staticmethod
    def plot1d(distribution, show=True, vmin=None, vmax=None, cmap=None):
        norm = mpl_colors.Normalize(
                vmin=vmin if vmin is not None else distribution.min(),
                vmax=vmax if vmin is not None else distribution.max()
        )
        cmap = mpl_cm.ScalarMappable(norm=norm, cmap=cmap or mpl_cm.get_cmap('jet'))
        cmap.set_array(distribution)
        c = [cmap.to_rgba(value) for value in distribution]  # defines the color

        fig, ax = plt.subplots()
        ax.scatter(np.arange(len(distribution)), distribution, c=c)
        fig.colorbar(cmap)
        if show: plt.show()
        return fig

    @staticmethod
    def plot2d(distribution, show=True):
        fig, ax = plt.subplots()
        img = ax.imshow(distribution, cmap='jet')
        fig.colorbar(img)
        if show: plt.show()
        return fig

    @staticmethod
    def plot3d(distribution, show=True):
        m, n, c = distribution.shape
        x, y, z = np.mgrid[:m, :n, :c]
        out = np.column_stack((x.ravel(), y.ravel(), z.ravel(), distribution.ravel()))
        x, y, z, values = np.array(list(zip(*out)))

        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')

        # Standalone colorbar, directly creating colorbar on fig results in strange artifacts.
        img = ax.scatter([0, 0], [0, 0], [0, 0], c=[0, 1], cmap=mpl_cm.get_cmap('jet'))
        img.set_visible = False
        fig.colorbar(img)

        ax.scatter(x, y, z, c=values, cmap=mpl_cm.get_cmap('jet'))
        if show: plt.show()
        return fig

Example

gaussian = Gaussian(size=(20,))
dist = gaussian.create(dim=1, centers=(1,)
gaussian.plot1d(dist, show=True)

Your problem

In order to get a solution that fits your question, the following code would work:

import numpy as np

arr = np.random.randint(0, 28, (25,))

gaussian = Gaussian(size=arr.shape)
centers = np.where(arr > 25.2)

distribution = gaussian.creates(dim=len(arr.shape), fwhm=2, centers=centers)
gaussian.plot(distribution, show=True)

For this the centers are determined by the condition arr > 25.2. Note that if there are no values, the next lines will crash. In order to get a width of 2 the value fwhm is put on 2, but you can alter this until you get the width that you want, or use Finding the full width half maximum of a peak.