3
$\begingroup$

I am dealing with the following function

fun[x_, z_] := a (1 - z^2) cos[x]^2 - (b z^2 + z - c).

My task is to generate the set of random pairs {x, z} such that the relation fun[x, z] > 0 holds. a, b and c are real numbers, -1 < z < 1 and 0 < x < 6. Any suggestion would be very much appreciated.

$\endgroup$
4
  • $\begingroup$ And do you want to choose these pairs uniformly in the relevant region? $\endgroup$
    – march
    Nov 22, 2015 at 5:11
  • $\begingroup$ No, it can be random either. $\endgroup$
    – Suro
    Nov 22, 2015 at 5:13
  • $\begingroup$ I'm sorry, I don't understand your comment. You could just choose a point uniformly (if that's the distribution from which you're pulling, which I assume) in the region specified be your inequality and check whether or not it satisfies fun[x,z] >0. If so, keep it. If not, throw it out. This is called rejection sampling (or the accept-reject Monte Carlo algorithm). $\endgroup$
    – march
    Nov 22, 2015 at 5:14
  • $\begingroup$ In case the point is uniformly, could you please help me to write an appropriate code $\endgroup$
    – Suro
    Nov 22, 2015 at 5:25

2 Answers 2

5
$\begingroup$

The function is given by

fun[x_, z_] := a (1 - z^2) Cos[x]^2 - (b z^2 + z - c)

We define a function that generates a random point in the x-z plane that satisfies fun[x, z] > 0:

accept[a0_, b0_, c0_] := 
  Block[{a = a0, b = b0, c = c0, x = RandomReal[{0, 6}], z = RandomReal[{-1, 1}]}
    , If[fun[x, z] > 0, {x, z}, accept[a, b, c]]
   ]

Important note: this definition is dangerous, because it is recursive. You should make sure that there are actually a good number of x-z pairs that satisfy fun[x, z] > 0 first before using this definition, because it runs the risk of running into the $RecursionLimit. As it is, it randomly generates a point, and if it satisfies the condition, it keeps it, and if it doesn't, then the function calls itself. Thus, it will eventually spit out a point if there is one to be found (or rather, if acceptable region is not too small).

Alternatively, you can define this function with a so-called "vanishing function" (named Nothing in Mathematica 10.1 and later) as

acceptReject[a0_, b0_, c0_] := 
  Block[{a = a0, b = b0, c = c0, x = RandomReal[{0, 6}], z = RandomReal[{-1, 1}]}
    , If[fun[x, z] > 0, {x, z}, ##&[]]
   ]

This will return literally nothing if the point doesn't satisfy the condition, and so when you make a Table of points, you will get only a list of points that satisfy the condition. However, the list of points won't necessarily be as many points as you specify, because it doesn't try to find another point if the random point that it chose doesn't fit the criterion; rather, it just moves on.

That said, we can generate a table of points via

points = Table[accept[3, -2, -2], {500}];

and we can plot the function along with these points to make sure that we are generating a uniform sample from the region:

Show[
 Plot3D[Evaluate[fun[x, z] /. {a -> 3, b -> -2, c -> -2}]
  , {x, 0, 6}, {z, -1, 1}
  , RegionFunction -> Function[#3 > 0]
  , Mesh -> None
  , PlotStyle -> Opacity[0.3]
  ]
 , Graphics3D[Point[{##, 0} & @@@ points]]
 ]

enter image description here

$\endgroup$
0
4
$\begingroup$
fun[a_, b_, c_, x_, z_] := a (1 - z^2) Cos[x]^2 - (b z^2 + z - c)
reg[a_, b_, c_] := 
 ImplicitRegion[
  fun[a, b, c, x, z] > 0 && -1 < z < 1 && 0 < x < 6, {x, z}]

Visualizations (using triple {3,-2,2} used in another answer):

rp = RegionPlot[reg[3, -2, -2]];
cp = ContourPlot[fun[3, -2, -2, x, z], {x, 0, 6}, {z, -1, 1}, 
   Contours -> {{0}}];
p3 = Plot3D[fun[3, -2, -2, x, z], {x, 0, 6}, {z, -1, 1}, 
   MeshFunctions -> {#3 &}, Mesh -> {{0}}];
Row[{rp,cp,p3}]

enter image description here

You can generate a list of random points for a given {a,b,c} using RandomPoint for region:

ListPlot[RandomPoint[reg[3, -2, -2], 1000, {{0, 6}, {-1, 1}}]]

enter image description here

$\endgroup$
2
  • $\begingroup$ I am thankful for you code. $\endgroup$
    – Suro
    Nov 23, 2015 at 23:49
  • $\begingroup$ Hello ubpdqn. Let me thank you again for your nice code. I tried to plot the reg function defined above using PolarPlot function but with no luck. Could you please help me in this? Thank you in advance. $\endgroup$
    – Suro
    Nov 26, 2015 at 20:20

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.