13

Why do I get a gap in the plot below and how can I fix it? (If you are interested in it, you can see a new related question: How to plot an implicit value funtion, which is also a little chanlenging)

Code:

Plot[InverseFunction[
   1 - 0.6*CDF[NormalDistribution[1, 0.3], #] - 
     0.4*CDF[NormalDistribution[3, 0.3], #] &][x], {x, 0, 1}, 
 PlotRange -> All, Exclusions -> None]

Result

Community
  • 1
ben
  • 173
  • 7

5 Answers5

9

The function you wish to plot happens to be the InverseSurvivalFunction of a MixtureDistribution with component distributions NormalDistribution[1, 0.3] and NormalDistribution[3, 0.3], and weights .6 and .4, respectively.

Using the built-in functions MixtureDistribution and InverseSurvivalFunction we get the desired result without an issue:

dist = MixtureDistribution[{6, 4}, {NormalDistribution[1, 0.3], NormalDistribution[3, 0.3]}]; 
Plot[InverseSurvivalFunction[dist, x], {x, 0, 1}]

enter image description here

You can also use InverseCDF to get the same output:

Plot[InverseCDF[dist, 1 - x], {x, 0, 1}]
(* same picture *)

Update: Addressing the question in the comments:

I do need to characterize D[x*InverseSurvivalFunction[dist, x], x]

Using the product rule and the inverse function theorem, you can define

derivative[x_] := (InverseSurvivalFunction[dist, y] +
     (x/(D[1 - CDF[dist, y], y] /. y -> InverseSurvivalFunction[dist, y]))) /. y -> x;

Column[Plot[#, {x, 0, 1}, ImageSize -> 400] & /@
  {x InverseSurvivalFunction[dist, x], Evaluate@derivative[x]}]

enter image description here

kglr
  • 394,356
  • 18
  • 477
  • 896
  • It's very helpful. Btw, can we do operations on InverseSurvialFunction, like D[x*InverseSurvivalFunction[dist, x], x]? It seems it doesn't work. Thanks! – ben Jan 02 '15 at 22:24
  • @ben, for the mixture distribution dist, D[x*InverseSurvivalFunction[dist, x], x] and D[x*InverseCDF[dist, 1-x], x] do not work; but they do work for "simpler" distributions: for example, D[InverseCDF[NormalDistribution[0, 1], x], x] and D[InverseSurvivalFunction[NormalDistribution[0, 1], x], x]. – kglr Jan 02 '15 at 23:07
  • I do need to characterize D[x*InverseSurvivalFunction[dist, x], x]. Do you have any suggestions which can remove the gap and also keep doing the derivation? Thanks. – ben Jan 02 '15 at 23:14
  • @ben, since we do get a closed form solution for D[1 - CDF[dist, x], x], maybe you can use Inverse Function Theorem? For example, (1/(D[1 - CDF[dist, x], x] /. x -> InverseSurvivalFunction[dist, x])) /. x -> .5 gives -2.00119. – kglr Jan 02 '15 at 23:50
  • @ben, please see the update... – kglr Jan 03 '15 at 00:32
  • that's really helpful. You are the expert in Mathematica^-^. One more question, given profit function yInverseSurvivalFunction[dist, x+y], for any x, I want to get the optimal value y(x) to maximize the profit function. Do you know how to plot y*(x)? Really appreciate your work and help. – ben Jan 03 '15 at 01:21
  • Let us continue this discussion in chat. – ben Jan 03 '15 at 01:23
7

You have complex numbers as a result in these ranges.

Check this:

Table[{i,
  InverseFunction[
    1 - 0.6*CDF[NormalDistribution[1, 0.3], #] -
      0.4*CDF[NormalDistribution[3, 0.3], #] &][i]}, {i, 0, 1, 0.01}]

The problem seems to be generated internally because of the real number in the function and also because of the fact that the plot uses a real number when sampling points for the plot. To see this behavior, look at these evolutions:

N[InverseFunction[
   1 - 6/10*CDF[NormalDistribution[1, 3/10], #] -
     4/10*CDF[NormalDistribution[3, 3/10], #] &][3/10]]

(*2.79765*)

InverseFunction[
  1 - 6/10*CDF[NormalDistribution[1, 3/10], #] -
    4/10*CDF[NormalDistribution[3, 3/10], #] &][0.3]

(*1.47655 + 0.475155 I*)
Sektor
  • 3,320
  • 7
  • 27
  • 36
Basheer Algohi
  • 19,917
  • 1
  • 31
  • 78
7

We get a little insight to the "bug" by writing the CDF in terms of Erfc:

 InverseFunction[
       1 - 0.6*CDF[NormalDistribution[1, 0.3], #] - 
           0.4*CDF[NormalDistribution[3, 0.3], #] &][.3]
 InverseFunction[
       1 - 0.6 1/2 Erfc[2.3570226039551585` (1 - #)] - 
           0.4 1/2 Erfc[2.3570226039551585` (3 - #)] &][.3]

1.47655 + 0.475155 I

1.47655 + 0.475155 I

(1 - 0.6*1/2 Erfc[2.3570226039551585` (1 - #)] - 
   0.4*1/2 Erfc[2.3570226039551585` (3 - #)]) &@%

0.3 - 4.85056*10^-16 I

we see that while CDF can not take an imaginary argument, Erfc can and the erroneous result is indeed a complex valued inverse of the function.

Edit -- a fix

it turns out we can use ConditionalExpression to force selection of a real inverse:

 Plot[ InverseFunction[
         ConditionalExpression[ 
           1 - 0.6*CDF[NormalDistribution[1, 0.3], #] - 
               0.4*CDF[NormalDistribution[3, 0.3], #], Element[#, Reals] ] &] @
              x , {x, 0, 1}]

enter image description here

george2079
  • 38,913
  • 1
  • 43
  • 110
  • Thanks a lot. Btw, can we do operations it, like Plot[D[xInverseFunction[ ConditionalExpression[ 1 - 0.6CDF[NormalDistribution[1, 0.3], #] - 0.4*CDF[NormalDistribution[3, 0.3], #], Element[#, Reals]] &]@ x, x], {x, 0, 1}] It seems it doesn't work. Thanks! – ben Jan 02 '15 at 22:42
6
f[x_] = 1 -
   0.6*CDF[NormalDistribution[1, 0.3], x] -
   0.4*CDF[NormalDistribution[3, 0.3], x];

ParametricPlot[{f[x], x}, {x, -1, 3},
 AspectRatio -> 1/GoldenRatio]

enter image description here

Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
  • Hi Bob, Thanks for your help. I will frequently use this function in the later analysis and same problem occurs in more complex functions. – ben Jan 02 '15 at 17:26
5

To start with, this should not be a problem for Mathematica. The non-inverse function is reasonably well-behaved:

f[u_] := 1-0.6 CDF[NormalDistribution[1, 0.3], u] - 0.4 CDF[NormalDistribution[3, 0.3], u];
Plot[f[u], {u, -1, 5}]

Mathematica graphics

(and one can verify, plotting or otherwise, that f' remains $<0$ so that f is 1-to-1). All I can suggest is to increase PlotPoints:

g = InverseFunction[f];
Plot[g[x], {x, 0, 1}, PlotRange -> {0, 4}, PlotPoints -> 30000]

nearly fills the gap around .3 but does nothing below .2, and takes over 400 sec here :(.

Mathematica graphics

You may want to report it to Wolfram.

A.G.
  • 4,362
  • 13
  • 18