PoissonGaussianLikelihoodWrapper.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-Gaussian distribution.
 * <p>
 * For each observed value the log-likelihood is computed from the Poisson-Gaussian distribution (a Poisson convolved
 * with a Gaussian). The mean of the Poisson distribution is set using the expected value generated by the provided
 * function. The 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.
 * The sum uses a non-normalised Poisson-Gaussian distribution for speed (see
 * {@link PoissonGaussianFunction#pseudoLikelihood(double)} ).
 */
public class PoissonGaussianLikelihoodWrapper extends LikelihoodWrapper
{
	private PoissonGaussianFunction p0, p1;
	final private boolean usePicard = false;

	/**
	 * 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
	 * @param a
	 *            The initial parameters for the function
	 * @param k
	 *            The observed values
	 * @param n
	 *            The number of observed values
	 * @param alpha
	 *            Inverse gain of the EMCCD chip
	 * @param s
	 *            The Gaussian standard deviation
	 */
	public PoissonGaussianLikelihoodWrapper(NonLinearFunction f, double[] a, double[] k, int n, double alpha, double s)
	{
		super(f, a, k, n);
		// We have two functions: One for no poisson counts and one for poisson counts.
		// They differ in their normalisation.
		p0 = PoissonGaussianFunction.createWithStandardDeviation(alpha, 0, s);
		p0.setUsePicardApproximation(usePicard);
		p1 = PoissonGaussianFunction.createWithStandardDeviation(alpha, 1, s);
		p1.setUsePicardApproximation(usePicard);
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.function.LikelihoodWrapper#computeLikelihood()
	 */
	public double computeLikelihood()
	{
		// Compute the negative log-likelihood to be minimised
		double ll = 0;
		for (int i = 0; i < n; i++)
		{
			final double e = f.eval(i);
			if (e <= 0)
				// Use the function with the non-Poisson normalisation
				ll -= p0.logProbability(data[i], e);
			else
				ll -= p1.logProbability(data[i], e);
		}
		return ll;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see gdsc.smlm.function.LikelihoodWrapper#computeLikelihood(int)
	 */
	public double computeLikelihood(int i)
	{
		final double e = f.eval(i);
		if (e <= 0)
			return p0.logProbability(data[i], e);
		return p1.logProbability(data[i], e);
	}

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