1
$\begingroup$

I would like to plot the solutions of an equation, for different values of a parameter. This is my code

cdf[x_] := CDF[NormalDistribution[0, 1], x];
pdf[x_] := cdf'[x];
cdf2[x_] = cdf[x]^20;
pdf2[x_] = cdf2'[x];

mix[x_, h_, p_] = h pdf[x - p] + (1 - h) pdf2[x - p];
ratio[x_, h_] = mix[x, h, 1]/mix[x, h, 0];
sol = NSolve[ratio[x, h] == 0.5 && (-5 < x < 5), x]
Plot[x /. sol, {h, 0, 1}]

For $h$ between 0 and 0.6 there should be three solutions, for $h$ above 0.6 just one. The plot, which runs in 5 minutes, returns only the solution for $h$ bigger than 0.6, while I believe I should get 3 different lines.

Can you help me fix this?

Thank you

$\endgroup$
1
$\begingroup$
cdf[x_] := CDF[NormalDistribution[], x];
pdf[x_] := cdf'[x];
cdf2[x_] = cdf[x]^20;
pdf2[x_] = cdf2'[x];

mix[x_, h_, p_] = h pdf[x - p] + (1 - h) pdf2[x - p] // Simplify;

ratio[x_, h_] = mix[x, h, 1]/mix[x, h, 0] // Simplify;

Clear[sol]

sol[h_?NumericQ] := NSolve[ratio[x, h] == 1/2 && (-5 < x < 5), x]

Generate data for a ListLinePlot. This is slow

hValues = 
  Join[Range[0, 0.1, 0.01], Range[0.15, 0.55, 0.05], 
   Range[0.551, 0.6, 0.001], Range[0.6, 1, 0.1]];

data = Table[{h, #} & /@ (x /. sol[h]), {h, hValues}];

data2 = GatherBy[data, Length];

ListLinePlot[{
  Rest@Flatten[data2[[1]], 1],
  Sequence @@ Transpose[data2[[2]]]},
 PlotLegends -> Placed[Automatic, {0.75, 0.45}]]

enter image description here

EDIT: If you want to plot for various powers, parametrize the functions with the power. Further, as stated in the comments below, "the better solution is to use ContourPlot as suggested by @UlrichNeumann."

Clear["Global`*"]

cdf[x_, pwr : _Integer?Positive : 1] := CDF[NormalDistribution[], x]^pwr;

pdf[x_, pwr : _Integer?Positive : 1] := (D[cdf[z, pwr], z] /. z -> x);

mix[x_, h_, p_, pwr : _Integer?Positive : 1] := 
  h pdf[x - p] + (1 - h) pdf[x - p, pwr];

ratio[x_, h_, pwr : _Integer?Positive : 1] :=
  mix[x, h, 1, pwr]/mix[x, h, 0, pwr];

ContourPlot[Evaluate[Simplify[ratio[x, h, #] == 1/2 & /@ {20, 10, 5}]],
 {h, 0, 1}, {x, -1/4, 9/4},
 FrameLabel -> (Style[#, 12, Bold] & /@ {"h", "x"}),
 ImageSize -> 324,
 MaxRecursion -> 5,
 PlotLegends ->
  Placed[StringForm["Power = ``", #] & /@ {20, 10, 5}, {0.75, 0.5}]]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for your answer, it's just great and very helpful. I've got just one question: why is there a little 'hole' between line 3 and 4? I would like to fix this if it were possible somehow. I tried to change $dh$ at .0005 but this results in an error. Thank you $\endgroup$ – Api Dec 16 '18 at 20:38
  • $\begingroup$ There is a gap because the gradient is very steep and the steps would need to be smaller as you suggest. The better solution is to use ContourPlot as suggested by @UlrichNeumann. $\endgroup$ – Bob Hanlon Dec 16 '18 at 21:16
  • $\begingroup$ Hey I was going through this code and I noted that when you change the power in $cdf2[x]$, it does not work anymore. I tried with using as exponent 5 and 10 and got the first time an error and the second time a wrong graph. Do you have any suggestions about how to fix it? Thank you $\endgroup$ – Api Jan 16 at 16:16
3
$\begingroup$

NSolve cannot solve the equation because the equation isn't numeric(depends on h)

For a first insight of the solution use

 ContourPlot[ratio[x, h] == 0.5, {x, -1, 3}, {h, 0, 1},MaxRecursion -> 5, FrameLabel -> {x, h}]

enter image description here

Try

sol[h_?NumericQ] := Values[NSolve[ratio[x, h] == 0.5 && (-5 < x < 5), x]//Flatten]

to get parameter dependent solutions.

sol[.5]
(*{-0.193146, 1.3878, 1.94506}*)
sol[.75]
(*{-0.193147}*)

Addenum

Unfortunately the solution cannot be plotted, because the number of solutions varies ...

Looking at the Contourplot it is possible to evaluate the Contourline using NDSolve:

First we need one point of the contourline, for example point h==0:

NSolve[{ratio[x, 0] == 1/2, -3 < x < 3}, x, Reals] (*{x -> 2.12759}*)

contourline doesn't change

H = NDSolveValue[{D[ratio[x, h[x]], x] == 0 ,h[2.127591638090098`] == 0}, h,{x, -1, 3}]

The left boundary of the solution range

H["Domain"][[1]] (* {-0.193105, 2.23873} *)
x0=%[[1]]  

fullfills ratio[x0,h]==1/2 x==x0 is also a solution .

Show[{Plot[H[x], {x, -1, 3}, PlotRange -> {0, 1}],ParametricPlot[{x0, h}, {h, 0, 1}]}, AxesLabel -> {x, h[x]}]

enter image description here

That's it, hope I could help!

$\endgroup$
  • $\begingroup$ I tried, but then this error was reported: SetDelayed::write: Tag List in {{x->-0.193133},{x->0.774654},{x->2.11263}}[h_?NumericQ] is Protected. $\endgroup$ – Api Dec 16 '18 at 18:20
  • $\begingroup$ Would you mind to produce the code for plotting it? I am sorry, I am an absolute beginner, I tried some guesses but I always get an empty graph $\endgroup$ – Api Dec 16 '18 at 18:26
  • 1
    $\begingroup$ The code is Plot[Evaluate[ sol[h]], {h, 0, 1}] but evaluation doesn't finish... $\endgroup$ – Ulrich Neumann Dec 16 '18 at 18:28
  • $\begingroup$ It takes 5 minutes and then it returns a graph which again shows just one of the solutions $\endgroup$ – Api Dec 16 '18 at 18:34
  • $\begingroup$ Sorry, I don't know why the Plot doesn't work as expected. $\endgroup$ – Ulrich Neumann Dec 16 '18 at 18:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.