# Plotting Quantile function of a mixture of lognormal distributions and its integral

I am trying to create a Mathematica function that is the Quantile function of a mixture distribution in which the constituent distributions are Lognormal. I want it to work inside a Manipulate where one is changing parameters of the mixture and the constituent distributions. I need to plot the quantile function and an integral of the quantile function inside the Manipulate.

If one tries to obtain a symbolic expression, Mathematica punts.

 Quantile[MixtureDistribution[{w,
1 - w}, {LogNormalDistribution[ma, sa],
LogNormalDistribution[mb, sb]}], q]


And I've tried inverting the SurvivalFunction of the mixture, but Mathematica complains that it doesn't know how to do that problem. (Probably because there actually is no closed form solution, which is precisely why it punted on the Quantile).

So, in order to approximate the quantile function, what I have ended up doing each time I change a parameter (such as w, ma, sa, mb or sb) is compute pairs of the survival function for values of q between 0 and 1 and q and then do an interpolation (order 1). It works but, since I have to do a bunch of these inside the Manipulate, it is about 1-2 orders of magnitude slower than I would like. Does anyone have insights on how I might do the computation faster? (Yes, the constituent functions really need to be lognormal)

It's also been suggested in a response to a prior draft of this question, that I unleash Mathematica's ability to compute the Quantile function numerically. When I try that method, it also works but is even slower. Here's an example of what I am trying to accomplish.

(*reparameterize lognormal dist*)
LogNormalDistribution3[mean_, median_] :=
LogNormalDistribution[Log[median],
Sqrt Sqrt[Log[mean] - Log[median]]];
(*manipulate*)
Manipulate[
Module[{m, fl, f, gl, g},
m = MixtureDistribution[{w, 1 - w}, {LogNormalDistribution3[ma, sa],
LogNormalDistribution3[mb, sb]}];
List[Range[0.0001 (*avoid infinity problem*), 1, 0.02],
Quantile[m, 1 - Range[0.0001, 1, 0.02]]]];
f = Interpolation[fl, InterpolationOrder -> 1] (* marginal cost*);
List[Range[0.02, 1, 0.02],
Table[Quiet@Integrate[f[x], {x, 0.0001, z}]/z, {z, 0.02, 1,
0.02}]]];
g=Interpolation[gl,InterpolationOrder->1];
(* actually need FindRoot code here *)
ListLogPlot[{fl, gl}, Joined -> True]],
{{w, 0.5}, 0, 1,
Appearance -> "Labeled"},
{{ma, 5000}, 2000, 8000,
Appearance -> "Labeled"},
{{sa, 1000}, 500, 1500,
Appearance -> "Labeled"},
{{mb, 5000}, 2000, 8000,
Appearance -> "Labeled"},
{{sb, 1000}, 500, 1500,
Appearance -> "Labeled"},
SynchronousUpdating -> True]


PS I guess I could create a mega interpolate in which I take a grid of al; plausible values of w, ma, sa, mb, sb and q and precompute, but in order to get the requisite precision over a large domain of parameter values, I think this would need to be a pretty huge object in memory.

Help!

• Oh, and to make matters worse, I end up having to use FunctionInterpolation to compute an integral of the quantile function over its domain. If anyone is interested why on earth I am doing this, one can characterize the marginal cost of insuring a pool as the quantile function of the distribution of costs of the pool. (Einav). The average cost for each quantile is the integral of the marginal cost function divided by the quantile value. – Seth Chandler Aug 24 '16 at 20:14
• Hi Seth. If you are trying to calculate something that is numerically intensive, then Manipulate doesn't seem like the appropriate medium / tool / presentational form to adopt, because it calculates in real time. Why not prepare the calculation overnight (or however long it takes), and then present it as an Animate of a collection of images. Always works best that way, for computationally intensive operations. – wolfies Aug 25 '16 at 19:05
• You may well be right, but I really like the interactive exploration format created by Manipulate. I have been hoping that I can use some cleverness to get these computations down so that they do work in real time, but it may be that what I am trying to do here is just too complicated. – Seth Chandler Aug 26 '16 at 13:12

I post this in case it is helpful:

fun[ma_, mb_, sa_, sb_, w_, x_] :=
1/2 w Erfc[(ma - Log[x])/(Sqrt sa)] +
1/2 (1 - w) Erfc[(mb - Log[x])/(Sqrt sb)]
q[ma_, mb_, sa_, sb_, w_, v_] :=
Quiet[u /. First@FindRoot[fun[ma, mb, sa, sb, w, u] - v, {u, 0.5}]]
Manipulate[
Module[{sol = q[ma, mb, sa, sb, w, v]},
Column[{
Plot[fun[ma, mb, sa, sb, w, x], {x, 0, 100}, PlotRange -> {0, 1},
Frame -> True, GridLines -> {{sol}, {v}},
Epilog -> {Red, PointSize[0.02], Point[{sol, v}]},
ImageSize -> 300], Row[ {"Quantile[", v, "]=", sol}] }
]], {ma, 1, 10}, {{mb, 2}, 1, 10}, {sa, 1, 10}, {sb, 1, 10}, {w, 0,
1}, {v, 0, 1}]

• fun is CDF of distribution
• FindRoot is being used to determine quantile
• dummy values as I do not really understand the context I apologize for any errors or misunderstanding.

• Thanks for thinking about this problem. Basically you've plotted the survival function or CDF and let one compute custom values of the inverse (the quantile function) via FindRoot. Unfortunately, both visually and computationally I really need to plot the quantile function and compute something that at least emulates the quantile function well enough to be able to integrate over it. – Seth Chandler Aug 25 '16 at 15:54
• @SethChandler apologies. Thank you for clarifying. I am currently 'off the grid' for a little while with no computer. Look forward to seeing answer :) – ubpdqn Aug 25 '16 at 21:32

Is this what you want?

Manipulate[
Quantile[MixtureDistribution[{w, 1 - w}, {LogNormalDistribution[ma, sa],
LogNormalDistribution[mb, sb]}], α],
{{α, 0.95}, 0.001, 0.999, Appearance -> "Labeled"},
{{w, 0.5}, 0, 1, Appearance -> "Labeled"},
{{ma, 1}, 0.1, 20, Appearance -> "Labeled"},
{{sa, 0.1}, 0.01, 2, Appearance -> "Labeled"},
{{mb, 1}, 0.1, 20, Appearance -> "Labeled"},
{{sb, 0.1}, 0.1, 2, Appearance -> "Labeled"}] • No. What I need is more like this ...Manipulate[ Plot[Quantile[ MixtureDistribution[{w, 1 - w}, {LogNormalDistribution[ma, sa], LogNormalDistribution[mb, sb]}], [Alpha]], {[Alpha], 0, 1}], {{[Alpha], 0.95}, 0.001, 0.999, Appearance -> "Labeled"}, {{w, 0.5}, 0, 1, Appearance -> "Labeled"}, {{ma, 1}, 0.1, 20, Appearance -> "Labeled"}, {{sa, 0.1}, 0.01, 2, Appearance -> "Labeled"}, {{mb, 1}, 0.1, 20, Appearance -> "Labeled"}, {{sb, 0.1}, 0.1, 2, Appearance -> "Labeled"}] – Seth Chandler Aug 24 '16 at 21:28
• Did you mean to put in a link? – JimB Aug 24 '16 at 21:29
• Sorry, @JimBaldwin, I have revised the reply to better explain the issue. – Seth Chandler Aug 24 '16 at 21:33
• You might want to add the term "plot" or "figure" in your question to be explicit. – JimB Aug 24 '16 at 21:41