vonMises.java

package statalign.model.ext.plugins.structalign;

import statalign.base.Utils;
import cern.jet.math.Bessel;

public class vonMises{
	/**
	 * 
	 * @param kappa - concentration parameter
	 * @param mean - mean angle
	 * @param angle - angle to be evaluated
	 * @return log density of distribution at @param angle
	 */
	static double logDensity(double kappa, double mean, double angle){
		return kappa * Math.cos(angle - mean) - Math.log(2 * Math.PI * Bessel.i0(kappa));
	}
	
	static double density(double kappa, double mean, double angle){
		return Math.exp(kappa * Math.cos(angle - mean)) / (2 * Math.PI * Bessel.i0(kappa));
	}
	
	/**
	 * 
	 * @param kappa - concentration parameter
	 * @param mean - mean angle
	 * @return - angle simulated from vonMises distribution with given parameters
	 */
	static double simulate(double kappa, double mean){
		double vm = 0, U1, U2, U3, a, b, c, f, r, z;
		boolean cont;
		a = 1 + Math.pow(1 + 4 * Math.pow(kappa, 2), 0.5);
		b = (a - Math.pow(2 * a, 0.5)) / (2 * kappa);
		r = (1 + Math.pow(b, 2)) / (2 * b);
		cont = true;
		while (cont) {
			U1 = Utils.generator.nextDouble();
			z = Math.cos(Math.PI * U1);
			f = (1 + r * z) / (r + z);
			c = kappa * (r - f);
			U2 = Utils.generator.nextDouble();
			if (c * (2 - c) - U2 > 0) {
				U3 = Utils.generator.nextDouble();
				vm = Math.signum(U3 - 0.5) * Math.acos(f) + mean;
				vm = vm % (2 * Math.PI);
				cont = false;
			}
			else {
				if (Math.log(c/U2) + 1 - c >= 0) {
					U3 = Utils.generator.nextDouble();
					vm = Math.signum(U3 - 0.5) * Math.acos(f) + mean;
					vm = vm % (2 * Math.PI);
					cont = false;
				}
			}
		}
		return vm;
	}
	// </vonMises>
}