PoissonGammaGaussianLikelihoodWrapper.java

package gdsc.smlm.function;

/*----------------------------------------------------------------------------- 
 * GDSC SMLM Software
 * 
 * Copyright (C) 2014 Alex Herbert
 * Genome Damage and Stability Centre
 * University of Sussex, UK
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *---------------------------------------------------------------------------*/

/**
 * This is a wrapper for any function to compute the negative log-likelihood assuming a Poisson-Gamma-Gaussian
 * distribution.
 * <p>
 * For each observed value the log-likelihood is computed from the Poisson-Gamma-Gaussian distribution (a Poisson
 * convolved with a Gamma distribution convolved with a Gaussian). The Poisson-Gamma distribution is derived
 * analytically in the paper Maximilian H Ulbrich & Ehud Y Isacoff, Nature Methods - 4, 319 - 321 (2007) to explain the
 * probability distribution of ADUs given a fixed photon level per pixel and set gain in an EM-CCD camera (The Poisson
 * distribution models the photon count and the Gamma distribution models the EM-gain). This is then numerically
 * convolved with a Gaussian distribution to model the read noise of the camera.
 * <p>
 * The distribution of Ulbrich & Isacoff has no analytical solution to the convolution with a Gaussian. However the
 * convolution with a Gaussian has the most effect when the counts are low. The Poisson-Gamma-Gaussian can be
 * approximated using the Poisson-Gamma and a partial convolution with a Gaussian at low counts. This method is provided
 * as Python source code within the supplementary information of the paper Mortensen, et al (2010) Nature Methods 7,
 * 377-383. This Java implementation is based on the Python code.
 * <P>
 * The mean of the Poisson distribution is set using the expected value generated by the provided function. The scale
 * (EM-gain) for the Gamma distribution and standard deviation of the Gaussian is fixed and set in the constructor. The
 * mean of the Gaussian is assumed to be zero.
 * <p>
 * The negative log-likelihood can be evaluated over the entire set of observed values or for a chosen observed value.
 */
public class PoissonGammaGaussianLikelihoodWrapper extends LikelihoodWrapper
{
	final private PoissonGammaGaussianFunction p;

	/**
	 * Initialise the function.
	 * <p>
	 * The input parameters must be the full parameters for the non-linear function. Only those parameters with gradient
	 * indices should be passed in to the functions to obtain the value (and gradient).
	 * 
	 * @param f
	 *            The function to be used to calculated the expected values (if modelling EMCCD data this should
	 *            evaluate the value without the bias)
	 * @param a
	 *            The initial parameters for the function
	 * @param k
	 *            The observed values (if using EMCCD data the bias should be subtracted)
	 * @param n
	 *            The number of observed values
	 * @param alpha
	 *            Inverse gain of the EMCCD chip
	 * @param s
	 *            The Gaussian standard deviation
	 */
	public PoissonGammaGaussianLikelihoodWrapper(NonLinearFunction f, double[] a, double[] k, int n, double alpha,
			double s)
	{
		super(f, a, k, n);
		p = new PoissonGammaGaussianFunction(alpha, s);
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.function.LikelihoodWrapper#computeLikelihood()
	 */
	public double computeLikelihood()
	{
		// Compute the negative log-likelihood
		double ll = 0;
		for (int i = 0; i < n; i++)
		{
			final double l = p.logLikelihood(data[i], f.eval(i));
			if (l == Double.NEGATIVE_INFINITY)
			{
				return Double.POSITIVE_INFINITY;
			}
			ll -= l;
		}
		return ll;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.function.LikelihoodWrapper#computeLikelihood(int)
	 */
	public double computeLikelihood(int i)
	{
		return -p.logLikelihood(data[i], f.eval(i));
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.function.LikelihoodWrapper#canComputeGradient()
	 */
	@Override
	public boolean canComputeGradient()
	{
		return false;
	}
}