Note: In the updated question, there is a genuine optimization problem, which is still amenable to the techniques of the original answer.
Original answer
I worked on the integral before I realized the optimization problem was trivial. Yet I thought I might as well share what I got.
It's sometimes hard to figure out when pursuing a symbolic solution will be faster and more productive than pursuing a numerical one. Numerical solutions are usually fast and approximate with sufficient accuracy, but in this case they are not fast. Normal distributions have been a subject of intense study, and so there's a chance that Mathematica can produce an exact solution. In this case that turns out to be worth trying.
First, when working with numerics in which you might have to adjust the precision (through WorkingPrecision
, say), it helps to use exact numbers for coefficients and parameters if available (or set the precision high enough with SetPrecision
). In this case the distributions to use are
BinormalDistribution[{0, 0}, {1, 1}, 1/2]
MultinormalDistribution[{0, 0, 0}, {{1, 1/2, 1/2}, {1/2, 1, 1/2}, {1/2, 1/2, 1}}]
Second, cap[y, z]
is the same as
CDF[BinormalDistribution[{0, 0}, {1, 1}, 1/2], -{y, z}]
and the built-in function is faster and more accurate.
Third, the triple integral can be solved exactly (with the indefinite limits y
and z
). The first (triple integral) took twice as long as the second (iterated integral):
Integrate[X*Joint2[X, Y, Z],
{Y, y, Infinity}, {Z, z, Infinity}, {X, -Infinity, Infinity}]
Fold[Integrate[#1, #2] &,
X * Joint2[X, Y, Z],
{{X, -Infinity, Infinity}, {Y, y, Infinity}, {Z, z, Infinity}}]
(* (E^(-(y^2/2)) (1 + Erf[(y - 2 z)/Sqrt[6]]) + E^(-(z^2/2)) Erfc[(2 y - z)/Sqrt[6]]) /
(4 Sqrt[2 Pi]) *)
(The timings were 30.770272
and 16.735236
resp.)
So a maximization problem can be set up with the integration having having been completed beforehand. Most likely it would be quick. I can't really put this to the test on the problem in the OP, since there is only one number in the domain that satisfies the constraint cap[y, 1] == 0.1
.
Update - Response to the additional problem
The new q
can still be done symbolically as above, which will save a lot of time over a triple NIntegrate
:
Fold[Integrate[#1, #2] &,
X*Joint2[X, Y, Z],
{{X, -Infinity, Infinity}, {Y, y, Infinity}, {Z, zstar[y], Infinity}}]
(*(E^(-(y^2/2)) (1+Erf[(y-2 zstar[y])/Sqrt[6]])
+E^(-(1/2) zstar[y]^2) \ Erfc[(2 y-zstar[y])/Sqrt[6]])/
(4 Sqrt[2 \[Pi]])*)
Note: it helps to do the X
integral first, probably in this case because it is a definite integral. This is true even for the triple Integrate
in the original answer above.
The next most expensive part of the calculation is FindRoot
. It can be sped up by about 50% by using InverseCDF
instead of CDF
:
z /. FindRoot[InverseCDF[Cond[z], 1 - .3] == y, {z, y - 0.1}]
But in this case FindRoot
can be replaced by Solve
, since Cond[z] == NormalDistribution[z/2, Sqrt[3]/2]
. Both of the following return the same thing (although the first generates an error message about using inverse functions:
z /. First@Solve[1 - CDF[Cond[z], y] == .3, z]
z /. First@Solve[InverseCDF[Cond[z], 1 - .3] == y, z] // Simplify
(* -0.908288 + 2. y *)
So zstar
can be defined by a linear equation (and could even be made exact, if desired).
With these changes NMaximize
takes a little time:
q[y_?NumberQ] := (E^(-(y^2/2)) (1 + Erf[(y - 2 zstar[y])/Sqrt[6]]) +
E^(-(1/2) zstar[y]^2) Erfc[(2 y - zstar[y])/Sqrt[6]]) / (4 Sqrt[2 Pi]);
zstar[y_?NumberQ] := Evaluate[z /. First@Solve[InverseCDF[Cond[z], 1 - .3] == y, z]];
NMaximize[{(1/2)*q[y1] + (1/2)*q[y2], cap[y1, y2] == .1}, {y1, y2}] // Timing
(* {2.683494, {0.0970647, {y1 -> 0.816697, y2 -> 0.816683}}} *)
cap[Y, Z] == CDF[BinormalDistribution[{0, 0}, {1, 1}, 1/2], {-Y, -Z}]
?CDF
is much faster thanNIntegrate
and probably more accurate. $\endgroup$cap[y, 1] == .1
has just one solution,y -> 0.462678
. So that will yield the maximum, as well the minimum :). $\endgroup$