# Conditional Expectation — How can Mathematica find a more general closed form?

Mathematica can't find a solution to this Expectation.

 Expectation[f \[Conditioned] 1 - (1 - f)^k >= p, {f, k} \[Distributed]
UniformDistribution[{{0, 1}, {2, 5}}],Assumptions->{p>0}]


But it can find solutions when I substitute specific rational numbers for p. Here's an example

 Expectation[
f \[Conditioned] 1 - (1 - f)^k >= 3/11, {f, k} \[Distributed]
UniformDistribution[{{0, 1}, {2, 5}}]]


It can even find solutions when I substitute transcendentals for p. Here's an example.

 Expectation[
f \[Conditioned] 1 - (1 - f)^k >= 1/\[Pi], {f, k} \[Distributed]
UniformDistribution[{{0, 1}, {2, 5}}]]


Is there some assumption I could add to the Expectation that would let Mathematica find a general solution? I've played around with various ideas such as p>0 but I can't find anything that works. Alternatively, can someone explain why Mathematica can't find a general solution?

P.S. If anyone cares, this question arises out of an attempt to understand subsidized insurance markets. Ultimately, I would like to use things fancier than UniformDistribution over which to compute the Expectation, but until I can get a close form solution to this simple case, more challenging cases look impractical.

-
I've done further work on this problem in light of Wolfie's response below. Thanks Wolfies for your efforts. They may well be correct. I could still use help from the community however. –  Seth Chandler Jul 27 '14 at 21:56
This does not really address closed form but some musings are here: ubpdqnmathematica.wordpress.com/2014/07/28/what-to-expect –  ubpdqn Jul 28 '14 at 5:44
Hi @SethChandler - I've added a second answer: the tweak you seek. :) –  wolfies Jul 28 '14 at 15:34

I figured out a way to answer a related simpler question, which is

 Expectation[
f \[Conditioned] f + \[Alpha]*k > p], {f, k} \[Distributed]
UniformDistribution[{{a, b}, {k1, k2}}]


Mathematica has trouble with this computation, but if you realize that it is the centroid of a polygon, you can use a feature of Mathematica to get a (long) answer

 RegionCentroid@
RegionIntersection[
ImplicitRegion[f + \[Alpha] k > p, {{f, a, b}, {k, k1, k2}}],
Rectangle[{a, k1}, {b, k2}]]


You then just take the first part of the answer, which provides the "f" or "x" coordinate.

I'm not sure if this is ugly or elegant and it is good enough to get me started on my project, but I would still enjoy an answer to the more challenging question posed above.

-

This is a very nice problem. If I may dispense with the notation for random variables as $f$ and $k$, and refer to them instead as $X$ and $Y$ ...

The Problem

Let $X \sim Uniform(0,1)$ and $Y \sim Uniform(2,5)$ be independent random variables, with joint pdf $f(x,y) = \frac13$:

   f = 1/3;    domain[f] = {{x, 0, 1}, {y, 2, 5}};


We seek a closed-form solution to:

$$E\big[\,X \, \, \big| \,\, \big\{(1 - X)^Y < 1-p \big\} \,\big] = \underset{\Omega}{{\int }}x \; \frac{f(x,y)}{P((x,y) \in \Omega)}$$

where $\Omega =\left\{(x,y): (1 - x)^y < 1-p)\right\}\,$ denotes the domain of support defined by the condition; integration is carried out over $\Omega$; and $0<p<1$.

To illustrate the dependency created by the condition, here is a plot of the domain of support $\Omega =\left\{(x,y): (1 - x)^y < 1-p)\right\}$, when $p=0.9$:

RegionPlot[(1-x)^y < .1, {x, 0, 1}, {y, 2, 5}, FrameLabel -> {x, y}]


The essence of the problem is to define the shaded space (the domain of support) in the above picture, given an arbitrary symbolic parameter $p$, in such a way that Mathematica can integrate over the domain of support $\Omega$, and so obtain a closed-form solution. Unfortunately, the usual simple trick of adding Boole[ (1-x)^y < 1-p ] into a manual integration does not work. In fact, even the simplest possible case, i.e.:

 Integrate[Boole[(1-x)^y < 1-p], {x, 0, 1}, {y, 2, 5}, Assumptions -> 0<p<1]


returns unevaluated input

So, if we require a closed-form symbolic solution for arbitrary parameter $p$, then other methods are needed.

Tackle the denominator first ...

The Denominator: $\quad P[(1-X)^Y < 1-p]$

Note that $P\big((1-X)^Y < q\big)$ where $q = 1-p$ is just the cdf of the distribution of $Z = (1-X)^Y$ ... so to get around the awkward non-rectangular domain of support, let us introduce the transformation:

$$Z = (1-X)^Y$$

Unfortunately, Mathematica 10 doesn't seem able to find this transformation by itself:

     CDF[TransformedDistribution[(1-x)^y,
{Distributed[x, UniformDistribution[{0, 1}]],
Distributed[y, UniformDistribution[{2, 5}]]}], q]


returns unevaluated input

No cigar. Let's try using the mathStatica package for Mathematica ...

Let $Z = (1-X)^Y$ and $W=Y$. Then, the joint pdf of $(Z,W)$ is:

Since $W=Y$, let us rather write this as the joint pdf of $(Z,Y)$, namely $g(z,y)$:

Then, the marginal pdf of $Z$ is say $h(z)$:

We seek $P\big((1-X)^Y < 1-p\big) = P(Z<1-p)$, which is:

Denominator: All done.

The Numerator:

We now wish to find the numerator, namely:

$$\underset{\Omega}{{\int }}x \; f(x,y) \quad \text{ where } \quad \Omega =\left\{(x,y): (1 - x)^y < 1-p)\right\}$$

Since $(X,Y)$ have joint pdf $f(x,y)$; and $(Z,Y)$ have joint pdf $g(z,y)$ (already derived); and since:

... the required integral is equivalent to finding:

$$\int _0^{1-p}\int _2^5 (1-z^{1/y}) \; g(z,y) \; dy \; dz$$

Note that we have converted a very awkward non-rectangular integration problem in $(X,Y)$ space into a neat, rectangular integral in $(Z,Y)$ space. Note too that the upper integral bound on $dz$ is $(1-p)$, since the imposed condition is $Z<1-p$. And so the numerator solution is:

Summary

The solution to $E\big[\,X \, \, \big| \,\, \big\{(1 - X)^Y < 1-p \big\} \,\big]\,$ is:

 sol = top/PP


... which yields the closed-form solution (copy and paste-able):

(2 + 10 (1 - p)^(1/5) - 5 (1 - p)^(2/5) - 4 Sqrt[1 - p] - 2 p + 2 Log[1-p] (-ExpIntegralEi[(1/5) Log[1-p]] + ExpIntegralEi[(2/5) Log[1-p]] + ExpIntegralEi[(1/2) Log[1-p]] - LogIntegral[1-p]))/ (2*(5*(1 - p)^(1/5) - 2 Sqrt[1-p] + (-ExpIntegralEi[(1/5) Log[1-p]] + ExpIntegralEi[(1/2) Log[1-p]]) Log[1-p]))

All done.

Here is a plot of the conditional expectation, as $p$ varies from 0 to 1:

Plot[sol, {p, 0, 1}, PlotRange -> {0, 1}, AxesLabel -> {p, "E[X | cond]"}]


When $p = 0$, there is no constraint, and so the conditional expectation is equal to the unconditional expectation $E[X] = \frac12$.

Finally, a quick numerical check: when $p = 0.9$:

 sol /. p -> 0.9


0.740363

which matches the numerical solution Mathematica can obtain, given a numerical value for $p$:

 Expectation[Conditioned[x, (1 - x)^y < 1 - 0.9],
Distributed[{x, y}, UniformDistribution[{{0, 1}, {2, 5}}]]]


0.740363

-
Can you see if MathStatica can solve this for the more general case Distributed[{x,y},UniformDistribution[{{a,b},{k1,k2}}]] instead of just Distributed[{x,y},UniformDistribution[{{0,1},{2,5}}]]. Also, do you happen to know if MathStatica plays with Mathematica 10? –  Seth Chandler Jul 28 '14 at 21:47
@SethChandler Hi Seth. Yes - mathStatica is fully compatible with v10. As to the general cases, have just tried: ...... f = 1/(d - c); domain[f] = {{x, 0, 1}, {y, c, d}} && {0 < c < d}; ... and that works perfectly, just following exactly the same steps. –  wolfies Jul 29 '14 at 16:38
@SethChandler If you also wish to generalise $X \sim Uniform(0,1)$ to $X \sim Uniform(a,b)$ ... since you have $1-X$, do you not need to assume that $0<x<1$ i.e. that $0\le a <b \le1$? –  wolfies Jul 29 '14 at 16:44

# Tweak Solution

Since Seth has provided 2 answers, I thought I might also put up another answer. My motivation for separating this from my original answer is that ...

• my original answer is a self-contained mathematical solution in transformations of random variables, essentially striving to side-step Mathematica's use of Boole which was not working, whereas ...

• the solution offered here is a simple Mathematica magic 'trick' that makes Boole work.

In my original answer, I noted that the problem involves calculations such as:

Integrate[Boole[(1-x)^y < 1-p], {x, 0, 1}, {y, 2, 5}, Assumptions -> 0<p<1]


returns unevaluated input

Seth has implored us to find an easier straightforward tweak solution.

Swap the order of integration

I just noticed that if one simply switches the order of integration:

qq = Integrate[Boole[(1-x)^y < 1-p], {y, 2, 5}, {x, 0, 1}, Assumptions -> 0<p<1]


5*(1-p)^(1/5) - 2*Sqrt[1-p] - ExpIntegralEi[(1/5)*Log[1-p]]*Log[1-p] + ExpIntegralEi[(1/2)*Log[1-p]]*Log[1-p]

... then Mathematica CAN suddenly magically DO the integral.

[ BTWqq/3 = $P\big((1-X)^Y < 1-p\big)$ = PP in my original answer. ]

That suggests that simply tricking Mma into switching the order of integration in Seth's original problem will make it too magically work, without any fuss or bother at all. Yum yums. And it does.

For example, by placing y \[Distributed] blah before x \[Distributed] bleh, so that x gets integrated out first:

 sol = Expectation[Conditioned[x, (1-x)^y < 1-p],
{y \[Distributed] UniformDistribution[{2, 5}],
x \[Distributed] UniformDistribution[{0, 1}]}]


... works:

... and which yields the same numerical test answer as before:

sol /. p -> .9


0.740363

-
Ingenious. But it appears to work only for the special case where the limits x[Distributed] UniformDistribution[{0,1}] and not for things like x[Distributed] UniformDistribution[{1/3,5/7}]. It is, however, more evidence that there is a bug in my solution and that reversing the order of integration shows some promise. Thanks for continuing to think about this problem (on which I am obsessing) –  Seth Chandler Jul 28 '14 at 16:50
@wolfies learned a lot from both your answers...interesting serendipitous observation but as Pasteur fortune favours the prepared mind –  ubpdqn Jul 29 '14 at 3:40

I spent some additional time on this problem and believe I MAY have solved it in a pretty general way with Mathematica. But I'm not sure and I still think there's more help needed on the problem. I'm going to use wolfies' substitution of the conventional x and y for my earlier used f and k.

Basically, what we are after is the expected value of x over a rectangular region subject to the constraint that x lies to the right of a curve:

 x > 1 - (1 - p)^y^(-1)


Alternatively, we can think of ourselves as trying to derive the expected value of x over that same rectangular region subject to the constraint that y lies above the curve

 Log[1 - p]/Log[1 - x]


Let's define the rectangular region as being bounded by a and b on the x-axis (0

We realize that, because x is uniformly distributed, the expected value of x for any value of y is simply the mean of two numbers: x1 and x2. So, if we could compute those two numbers as a function of y (and some other stuff), we could compute the mean and then convert this two-dimensional problem into a one-dimensional problem in which we compute ...

 Expectation[Mean[{f1[y],f2[y]}], Distributed[k,UniformDistribution[{k1, k2}]]]


The f2[y] part is simple: the right hand of the interval of interest is just b. The f1[y] part is just a little tougher. Either f1[y] is a or f1[y] is 1 - (1 - p)^y^(-1). It just depends on whether y>Log[1-p]/Log[1-a]. So the problem reduces to the following expectation:

 Expectation[
1/2 (b + If[k > Log[1 - p]/Log[1 - a], a, 1 - (1 - p)^(1/k)]),
k \[Distributed] UniformDistribution[{k1, k2}]]


The only other catch is that k1, the lower bound of the rectangular region has to be redefined to be the greater of your original k1 and Log[1-p]/Log[1-b]. So we now go get some coffee and ask Mathematica to evaluate the following problem

 FullSimplify[
Expectation[
1/2 (b + If[k > Log[1 - p]/Log[1 - a], a, 1 - (1 - p)^(1/k)]),
k \[Distributed] UniformDistribution[{k1, k2}]] /.
k1 -> Max[k1, Log[1 - p]/Log[1 - b]], {0 < a < 1, k2 > k1,
1 > b > a, k1 > 1, a < p < b}]


It turns out that Mathematica does know how to solve this problem. We get the following gory Piecewise expression:

 Piecewise[{{(a + b)/2,
Log[1 - a]*Max[k1, Log[1 - p]/Log[1 - b]] <= Log[1 - p] &&
k2*Log[1 - a] < Log[1 - p] &&
Max[k1, Log[1 - p]/Log[1 - b]] < k2},
{((a + b)*k2 +
Log[1 - p]*(Gamma[
0, -(Log[1 - p]/Max[k1, Log[1 - p]/Log[1 - b]])] +
LogIntegral[1 - a]) + (-1 -
b + (1 - p)^(1/Max[k1, Log[1 - p]/Log[1 - b]]))*
Max[k1, Log[1 - p]/Log[1 - b]])/(2*(k2 -
Max[k1, Log[1 - p]/Log[1 - b]])),
Max[k1, Log[1 - p]/Log[1 - b]] < k2 &&
Log[1 - a]*Max[k1, Log[1 - p]/Log[1 - b]] > Log[1 - p] &&
k2*Log[1 - a] < Log[1 - p]},
{(k2*(1 + b - (1 - p)^(1/k2)) + (ExpIntegralEi[Log[1 - p]/k2] -
ExpIntegralEi[Log[1 - p]/Max[k1, Log[1 - p]/Log[1 - b]]])*
Log[1 - p] + (-1 -
b + (1 - p)^(1/Max[k1, Log[1 - p]/Log[1 - b]]))*
Max[k1, Log[1 - p]/Log[1 - b]])/(2*(k2 -
Max[k1, Log[1 - p]/Log[1 - b]])),
Max[k1, Log[1 - p]/Log[1 - b]] < k2 &&
k2*Log[1 - a] >= Log[1 - p] &&
Log[1 - a]*Max[k1, Log[1 - p]/Log[1 - b]] >= Log[1 - p]},
{(k2*(1 + a +
2*b - (1 - p)^(1/k2)) + (ExpIntegralEi[Log[1 - p]/k2] -
ExpIntegralEi[Log[1 - p]/Max[k1, Log[1 - p]/Log[1 - b]]])*
Log[1 - p] + (-1 - a -
2*b + (1 - p)^(1/Max[k1, Log[1 - p]/Log[1 - b]]))*
Max[k1, Log[1 - p]/Log[1 - b]])/(2*(k2 -
Max[k1, Log[1 - p]/Log[1 - b]])),
Max[k1, Log[1 - p]/Log[1 - b]] < k2 &&
k2*Log[1 - a] >= Log[1 - p] &&
Log[1 - a]*Max[k1, Log[1 - p]/Log[1 - b]] < Log[1 - p]}}, 0]


I wanted to check whether this answer (soln) was consistent with Wolfies, which was derived with the help of the MathStatica package. So, I wrote

 Refine[soln /. {a -> 0, b -> 1, k1 -> 2, k2 -> 5}, 0 < p < 1]


And got ...

 1/6 (5 (2 - (1 - p)^(1/5)) +
2 (-2 + Sqrt[1 - p]) + (ExpIntegralEi[1/5 Log[1 - p]] -
ExpIntegralEi[1/2 Log[1 - p]]) Log[1 - p])


This is similar to Wolfies' (it has those wierdo ExpIntegralEi things) but NOT identical. And when I evaluate his solution and mine for a=0,b=1,k1=2,k2=5,p=0.9 I get 0.748267, which is really close to his 0.740363 but NOT identical.

So, I am still left with a couple of mysteries:

(a) am I doing this correctly? (b) why is Wolfies' solution slightly different than mine? and (c) why can't Mathematica do this problem more directly without requiring heroics (a/k/a thinking) on the part of the user. It would be nice if Mathematica could actually do the following:

 Expectation[
x \[Conditioned] x >= 1 - (1 - p)^y^(-1), {x, y} \[Distributed]
UniformDistribution[{{a, b}, {k1, k2}}],
Assumptions -> {0 < a < b < 1, 1 < k1 < k2, 0 < p < 1}]


More help still appreciated.

-
(a) not quite right ; (b) because his is right; (c) I don't know. I'm still working on this problem. –  Seth Chandler Jul 28 '14 at 22:17