java - 计算最小化另一个函数的正常函数的反函数

Question

我实现了一个方法 (Normal.compute())，用于计算两个法线函数的总和:

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));
    }

}

现在固定了一个值z和一个点(x1,y1)，我将检索最近的(就距离而言)点，其值 >Normal.compute() = z。其中最近的位置由 distance() 计算。

所以我需要的是计算反函数并最小化距离，但我不知道如何以编程方式实现。

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 }
}

我尝试过使用 apache common math 或 colt，但它们似乎没有帮助。

这不是练习，所以如果已经完成了某些操作，我可以使用库。

Answer 1

您可以将其表述为优化嵌套在查找根中的函数。

优化问题:找到max_t normal.compute(foo(t))，其中foo(t)参数化一个以r为中心半径的圆在(x1, y1)上，即foo(t) = [x1 + r cos(t), y1 + r sin(t)]。 (如果 normal.compute(x1, y1) 小于 z，则使用 max。否则使用 min。)

寻根问题:找到r，使得优化问题的解 = z。

我通过查看 normal.compute(x, y) 的等高线图并思考 (x1, y1) 周围的圆圈得到了这个想法。当您绕圈行走时，normal.compute 的值会上下变化。您希望圆的最高点或最低点恰好是 z。希望这会有所帮助。

编辑:我喜欢这个问题，所以我构建了一个解决方案。这是 Maxima [1] 计算机代数系统的脚本。请参阅代码中的注释。玩得开心。

/* 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 */