Calculate inverse of normal function that minimizes another function

Question

I have implemented a method (Normal.compute()) that calculates the sum of two normal function:

public class Normal {
    private double mu1, mu2;
    private double sigma1, sigma2;

    public double compute(double x, double y) {
        return normal(x,mu1,sigma1) + normal(y,mu2,sigma2);
    }

    public double normal(double x, double mu, double sigma) {
        return Math.pow(Math.E, (-1*Math.pow((x-mu),2)) /   (2* Math.pow(sigma,2))
                        /
                        sigma * Math.sqrt(2*Math.PI)
        );
    }

    public static double distance(double ax1, double ax2, double bx1, double bx2) {
        return Math.sqrt(Math.pow((bx1-ax1),2)+Math.pow((bx2-ax2),2));
    }

}

Now fixed a value z and a point (x1,y1), I would retrieve the nearest (in terms of distance) point for which the value of Normal.compute() = z. Where nearest is computed by distance().

So what I need is to calculate the inverse function and minimize the distance, but I don't know how do it programmatically.

public double[] inverseNearest(double z, double x1, double x2) {
    // K = Set of (x,y) such that compute(x,y) = z

    // return argmin { distance(xk, yk, x1, x2) for each (xk, yk) in K }
}

I've tried to use apache common math or colt, but they don't seem helpful.

This is not an exercise, so I can use a library if there is something already done.

Answer 1

You can phrase this as optimizing a function nested inside finding a root.

Optimization problem: find max_t normal.compute(foo(t)) where foo(t) parametrizes a circle of radius r centered on (x1, y1), i.e. foo(t) = [x1 + r cos(t), y1 + r sin(t)]. (Use max if normal.compute(x1, y1) is less than z. Otherwise use min.)

Root finding problem: find r such that the solution to the optimization problem = z.

I got this idea by looking at a contour plot of normal.compute(x, y) and thinking about circles around (x1, y1). The value of normal.compute goes up and down as you walk around the circle. You want the circle such that its highest or lowest point is exactly z. Hope this helps.

EDIT: I like this question, so I've constructed a solution. Here is a script for the Maxima [1] computer algebra system. See comments in the code. Have fun.

/* solve problem stated in:
 * http://stackoverflow.com/questions/22099321/calculate-inverse-of-normal-function-that-minimizes-another-function
 *
 * copyright 2014 by Robert Dodier
 * I release this work under terms of the GNU General Public License
 *
 * how to:
 *
 * (1) assign values to m1, s1, m2, s2, x1, y1, and z0
 * (2) batch("foo.mac");
 *
 * example:
 *
 * [m1, s1, m2, s2, x1, y1, z0] : [-2, 1.8, 1.3, 2.4, 1.82, -0.24, 0.3];
 * batch ("foo.mac");
 *  => [x0, y0] = [- .4249300563696112, 0.172672148095035]
 *
 * after that, try:
 *
 * set_plot_option ([same_xy, true]);
 * load (implicit_plot);
 * implicit_plot ([F (x, y) = z0, (x - x1)^2 + (y - y1)^2 = r0^2], [x, -5, 5], [y, -5, 5]), numer;
 *
 * should show z0 contour just touching the circle of radius r0 centered on [x1, y1]
 */

load (distrib);
ratprint : false;

/* m1, s1, m2, s2 must be assigned values */

F (x, y) := pdf_normal (x, m1, s1) + pdf_normal (y, m2, s2);

/* x1, y1 must be assigned values */

G (r, x1, y1) := fmax_circular (lambda ([t], F (x1 + r * cos (t), y1 + r * sin (t))));

fmax_circular (f) := lmax (map (f, fargmax_circular (f)));

fargmax_circular (f) := block ([n : 17, u0, u2],
  map (f, makelist (i * 2 * float (%pi) / n, i, 0, n)),
  ev (sublist_indices (%%, lambda ([u], u = u_max)), u_max = lmax (%%)),
  map (lambda ([i], [u0, u2] : [(i - 1) * 2 * float (%pi) / n, (i + 1) * 2 * float (%pi) / n], fargmax1 (f, u0, u2)), %%));

/* golden section search */

fargmax1 (f, u0, u2) := block ([tol : 1e-2, u1 : u0 + (u2 - u0) / float (%phi)],
  while u2 - u0 > tol
    do block ([x],
      if u2 - u1 > u1 - u0
        then x : u1 + (1 - 1/float (%phi)) * (u2 - u1)
        else x : u1 - (1 - 1/float (%phi)) * (u1 - u0),
      if f(x) > f(u1)
        then /* accept interval containing x */
          if u2 - u1 > u1 - u0
            then [u0, u1, u2] : [u1, x, u2]
            else [u0, u1, u2] : [u0, x, u1]
        else /* reject interval containing x */
          if u2 - u1 > u1 - u0
            then [u0, u1, u2] : [u0, u1, x]
            else [u0, u1, u2] : [x, u1, u2]),
  u0 + (u1 - u0) / 2);

/* z0 must be assigned a value */

r0 : find_root (lambda ([r], G (r, x1, y1) - z), r, 0.001, 5), z = z0, numer;

/* fargmax_circular returns a list -- assume it's just one element
 * not guaranteed to work -- [t0] : ... fails when rhs has 2 or more elements, oh well
 */
[t0] : fargmax_circular (lambda ([t], F (x1 + r0 * cos (t), y1 + r0 * sin (t))));

/* [x0, y0] is the point on z0 contour of F, nearest to [x1, y1]
 * r0 is distance from [x1, y1] to [x0, y0]
 */

[x0, y0] : [x1 + r0 * cos (t0), y1 + r0 * sin (t0)];

F (x0, y0), numer; /* should be equal to z0 */

Answer 2

normal(x,mu1,sigma1) + normal(y,mu2,sigma2)

That is not the correct way to compute the probability of the sum of two normally distributed random variables. If your x and y are drawn from two independent normal distributions, call them X and Y, then the sum X+Y is in turn normally distributed. The mean of Z=X+Y is the sum of the two means, and the square of the variance is the sum of the squares of the two variances.

In other words, you should be using

normal(x+y,mu1+mu2,Math.sqrt(Math.pow(sigma1,2)+Math.pow(sigma,2)))