Contour plots of a two-variable function depicting the values of its integral

Question

Given a two-variable function f[x,y], how to get the contour plots corresponding to several values of the integral of the original function. For example, I would like to plot the contours $C_j$ of the domains $S_j$, in the $xy$ plane, such that

$$ \int_{S_j} dxdy f(x,y) = \alpha_j $$

for $\alpha_j = \alpha_{j-1} + \epsilon$, with a fixed step $\epsilon$.

Additionally, I would like to be able to choose a value of $\alpha$ and plot the contour of the corresponding domain.

My difficulty is the following. I don't know how to make a condition inside ContourPlot such that it tells that the contours should correspond to the values of the integrals.

A possible application would be to determine the region corresponding to 68.3% of the events of a two-variable Gaussian distribution. This is called a 1-sigma region.

Answer 1

Here is a difficult running example:

f[x_, y_] := Abs[x + y]^(1/2) Cos[x^2 + y^2/2]^2;
Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, Boxed -> False]

I interpret the question in this way: given a value $z$ in the range of $f$, we need to be able to compute the integral of $f$ over the set of $R_z=\{(x,y)\ |\ f(x,y) \le z\}$. Here is a general way:

i[f_, z_, xlim_: {-Infinity, Infinity}, ylim_: {-Infinity, Infinity}] := 
  Block[{x, y, g = Function[{x, y}, u = f[x, y]; u Boole[u <= z]]}, 
  Chop[NIntegrate[g[x, y], Evaluate@Prepend[xlim, x], Evaluate@Prepend[ylim, y]]]];

Its arguments are $f$, $z$, and the endpoints of a rectangle containing the domain of $f$ (defaulting to all of $\mathbb{R}^2$. For example,

i[f, 0.5, {-2, 2}, {-2, 2}]

returns $1.69395$, which is the numerical estimate of the integral $\int_{R_{0.5}} f(x,y) dx dy$.

The question asks how to invert $i$. To this end, let's compute the values of $i(f,z)$ for a carefully chosen range of values of $z$ matching the characteristics of $f$ (as gleaned by inspecting the preceding plot):

x = Range[2/10, 13/10, 1/10]^1.5;
data = ParallelMap[i[f, #, {-2, 2}, {-2, 2}] &, x]

Because $f$ is a somewhat difficult function to integrate, this takes some time and raises a lot of warnings. You can improve the accuracy according to needs and time available. But let's continue. Here is what we just computed:

Show[ListPlot[{x, data}\[Transpose], Joined -> True, PlotStyle -> Thick], 
  ListPlot[{x, data}\[Transpose], PlotStyle -> PointSize[0.015]]]

That's smooth enough for this illustration; here's the interpolating function of its inverse:

h = Interpolation[{data, x}\[Transpose]]

With it we can quickly produce any of the desired contours. Here, we show them for $z=1,2,\ldots, 7$:

ContourPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, Contours -> h /@ Range[1, 7]]

The contours themselves are labeled (as tooltips, not visible in this figure) with the $h(z_i)$, giving the height cutoff for $f$ corresponding to each contour.

Obviously there are some errors near the extreme values of $z$, related to the finite accuracy with which the integrals have been computed. Recompute them with higher accuracy if desired. For $f$ with sufficiently simple expressions, NIntegrate can be replaced by Integrate in i to achieve perfect accuracy and the interpolation in h can even be replaced by a call to Solve.