12
$\begingroup$

I'm writing a program that i need 1000 real number between -3 and 6. the code that I used for generating and saving them was

l = RandomReal[{-3, 6}, 1000]

But I want to be sure none of these numbers are equal to 3.How can I do that in a short and nice way? How can I generate these numbers in a way that I don't need to check for 3 being in list?

$\endgroup$
  • 5
    $\begingroup$ It might help to explain why you wan't this. In theory the odds of generating 3 exactly is infinitesimally small (it may in fact be impossible depending on the specifics of the generator) and as such not really worth bothering about. If you are rounding off a some point this may be different. $\endgroup$ – user22695 Nov 28 '14 at 5:22
  • 1
    $\begingroup$ I might ask how sure you need to be. The probability of generating a machine number that is detected as equal to 3. is very very small. Somewhere on the order of 10^-12 when you only generate 1000 of them in that range. It is about two orders of magnitude smaller if you are checking for equivalence (===). $\endgroup$ – Andy Ross Nov 28 '14 at 5:23
11
$\begingroup$
Plot[PDF[MixtureDistribution[{6, 3},
        {UniformDistribution[{-3, 3 - $MachineEpsilon}], 
         UniformDistribution[{3 + $MachineEpsilon, 6}]}], x], {x, -3, 6}, 
         Filling -> Axis]

Mathematica graphics

To generate your data:

RandomVariate[ MixtureDistribution[{6,3}, 
              {UniformDistribution[{-3, 3 - $MachineEpsilon}], 
           UniformDistribution[{3 + $MachineEpsilon, 6}]}], 500]
$\endgroup$
  • 2
    $\begingroup$ I do not believe this works. Sorry, but -1. Explanation to follow. $\endgroup$ – Mr.Wizard Nov 28 '14 at 12:56
  • $\begingroup$ No idea whether it work ... but I do not believe it is necessary. Just using RandomReal should be just fine. $\endgroup$ – wolfies Nov 28 '14 at 15:19
  • $\begingroup$ Nice!I just wanna vote a +1. $\endgroup$ – yode Dec 27 '15 at 0:56
9
$\begingroup$

I do not believe that belisarius's answer, as written, is correct, at least in Mathematica 10.0.1 under Windows.

In a simplified example we can see that in machine precision a three is still generated despite subtracting $MachineEpsilon from 3 in the range. I narrow the range make this readily apparent:

ϵ = $MachineEpsilon;

SeedRandom[1]
First /@ RealDigits @ RandomReal[{3 - 5 ϵ, 3 - ϵ}, 5]
{{3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
 {2,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9},
 {3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
 {2,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9},
 {2,9,9,9,9,9,9,9,9,9,9,9,9,9,9,9}}

We can see that at least in Mathematica 10.0.1 3 - ϵ is not distinct from 3.:

{3., 3 - ϵ} // FullForm
Order[3., 3 - ϵ]
List[3.`, 3.`]

0

The longer formulation is not even accepted:

RandomVariate[
 MixtureDistribution[{6, 3},
  {UniformDistribution[{3 - 5 ϵ, 3 - ϵ}], 
   UniformDistribution[{3 + ϵ, 3 + 5 ϵ}]}], 500]
UniformDistribution::lss: Parameter 2.999999999999999` at
position {1,1} in UniformDistribution[{3.,3.}] is expected to be less
than 3.`. 

UniformDistribution::lss: Parameter 2.999999999999999` at position
{1,1} in UniformDistribution[{3.,3.}] is expected to be less than 3.`.

This is despite 3 - 5 ϵ and 3 - ϵ being distinct. If I expand the range sufficiently the input is accepted, but as one might expect following the behavior illustrated earlier three is included in the result:

SeedRandom[1]

RandomVariate[
  MixtureDistribution[{6, 3},
   {UniformDistribution[{3 - 500 ϵ, 3 - ϵ}], 
    UniformDistribution[{3 + ϵ, 3 + 500 ϵ}]}], 5] // RealDigits
{{{3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 4}, 1},
 {{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 8, 9, 6}, 1},
 {{3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7}, 1},
 {{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 4, 9}, 1},
 {{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 1, 4}, 1}}

For the simplified example it appears sufficient to use an offset of 2 ϵ as this is distinct from three:

SeedRandom[1]
RandomReal[{3 - 5 ϵ, 3 - 2 ϵ}, 500000] // RealDigits // Tally
{{{{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9}, 1}, 500000}}

However for the longer formulation this too fails. In fact even using an offset of 100 ϵ is insufficient!

SeedRandom[1]

RandomVariate[
  MixtureDistribution[{6, 3},
   {UniformDistribution[{3 - 500 ϵ, 3 - 100 ϵ}], 
    UniformDistribution[{3 + 100 ϵ, 3 + 500 ϵ}]}], 5] // RealDigits
{{{3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 7}, 1},
 {{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 8, 9, 4}, 1},
 {{3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 4}, 1},
 {{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 3, 7}, 1},
 {{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 0, 9}, 1}}

As a final example observe that I can recreate belisarius's plot without using any offset at all:

Plot[
 PDF[
  MixtureDistribution[{6, 3},
   {UniformDistribution[{-3, 3}], 
    UniformDistribution[{3, 6}]}], x],
 {x, -3, 6},
 Filling -> Axis
]

enter image description here


Further problems

Contrary to what is stated (and illustrated) above an offset of 2 ϵ may not be sufficient for RandomReal either. There is apparently instability in this process if not a bug. While continuing to experiment I got a different result:

ϵ = $MachineEpsilon;

SeedRandom[1]

RandomReal[{3 - 5 ϵ, 3 - 2 ϵ}, 500000] // RealDigits // Tally
{{{{3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1}, 124641},
 {{{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9}, 1}, 375359}}

Yet the same code in a fresh kernel again yields:

{{{{2, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9}, 1}, 500000}}

Proposed solution

Since I am unable to recommend any method to directly avoid generating a specific value within a range, especially in light of the apparently unstable behavior I am experiencing from RandomReal, I can only suggest a filtering approach which will be admittedly slow. Here with integers for an easier illustration:

RandomInteger[{1, 5}, 50] //. {4 :> RandomInteger[{1, 5}]}
{1, 2, 3, 3, 1, 5, 2, 5, 1, 5, 3, 2, 5, 2, 3, 3, 2, 1, 3, 2, 5, 5, 3, 2,
 5, 5, 1, 1, 1, 3, 2, 5, 1, 2, 3, 1, 2, 1, 5, 1, 2, 3, 1, 1, 5, 1, 3, 2, 3, 1}

This should robustly assure that no 4 appears in the output, and that the output is always of length 50, no matter how many 4's appear in the initial result. The same method can be applied to real values. This is a great deal of overhead for an unlikely event however, so consider your alternatives carefully.

$\endgroup$
  • 1
    $\begingroup$ I didn't check this or even look too hard at this question/answers. But, it should be 3*(1 - $MachineEpsilon) and 3*(1 + $MachineEpsilon), not as written in belisarius's answer. After that, it should work correctly (if there are no FP bugs in the distribution implementations, &c.). $\endgroup$ – Oleksandr R. Nov 28 '14 at 13:26
  • $\begingroup$ m = $MachineEpsilon; RandomVariate[ MixtureDistribution[{1, 1}, {UniformDistribution[{-2 m, -m}], UniformDistribution[{m, 2 m}]}], 500] $\endgroup$ – Dr. belisarius Nov 28 '14 at 13:32
  • 1
    $\begingroup$ @belisarius I don't follow? Anyway, reading this answer more carefully, I see that (i.e., the error shows that) Andy (or whoever else wrote or helped to write these) has been bitten by Internal`$EqualTolerance and friends. These are sufficiently obscure that it wouldn't surprise me at all that there is at least one person working at WRI who is writing numerical code and doesn't know what they do, or considers it an irrelevant corner case. $\endgroup$ – Oleksandr R. Nov 28 '14 at 13:43
  • 1
    $\begingroup$ Amusingly(?), SetPrecision[Union@RandomVariate[UniformDistribution[{1.-$MachineEpsilon,1.}],100],Infinity] shows that Mathematica also does not seem to be quite sure which definition of $\epsilon$ it wants to use--is it 1/2 ULP (LAPACK convention), or 1 ULP (ISO C convention)? $\endgroup$ – Oleksandr R. Nov 28 '14 at 14:05
  • 2
    $\begingroup$ @belisarius $MachineEpsilon is normalized to 1, i.e. it's a relative, not an absolute, difference. There do exist plenty of valid FP numbers smaller than it, but not any two that vary by a fraction less than it. $\endgroup$ – Oleksandr R. Nov 28 '14 at 14:07
5
$\begingroup$

Theoretically, for any continuous distribution ... which includes your Uniform distribution ... $P(X=3) = 0$. So, theoretically, your question leaves you with nothing to worry about.

Practically, for RandomReal to hit a perfect 3 ... should be equally impossible. Perhaps, as a matter of machine-precision, you are worried about whether mma might hypothetically assert a real is equal to an integer. It should not.

But, if it really, really bothers you, ... simply generate your pseudo-random data ... check if it contains a 3 ... and if it does, delete it.

Just out of interest, here are 100 million pseudo-random numbers drawn between 1.9999 and 2.00001 ...

data = RandomReal[{1.9999, 2.00001}, 10^8];

and there is not one integer 2:

MemberQ[data, 2]

False

Nor should there be.

$\endgroup$
  • $\begingroup$ But perhaps there is a 2.0. Yes, the real number 2 is actually thought of by Mathematica as a uniform distribution having width 3 $MachineEpsilon, and yes, "equality" is not really defined for distributions, but this is a unique oddity (some would say perversion) of Mathematica and probably not what the OP is asking about. Practically speaking, the calculation is done over a finite set of rationals, not reals or distributions or anything else. $\endgroup$ – Oleksandr R. Nov 28 '14 at 15:53
  • 2
    $\begingroup$ Yes. P(X=3)=0. This is the correct answer to this question. $\endgroup$ – dnet Nov 28 '14 at 17:03
  • 1
    $\begingroup$ The OP clearly states that they do not want numbers "equal to 3" in the result. Integers are a subset of the reals and 3 is equal to 3.0 (but it is not a pattern match in Mathematica, as you surely know). You are obviously correct if Mathematica really draws from a uniform distribution over the reals, but given that (like every other computer system) it does this only using the approximation of floating point arithmetic, then either your answer is wishful thinking, or your math is faulty. $\endgroup$ – Oleksandr R. Nov 28 '14 at 18:08
  • 1
    $\begingroup$ data = RandomReal[{1.99999999, 2.00000001}, 10^8]; MemberQ[data, 2.] (* -> True *) $\endgroup$ – Oleksandr R. Nov 28 '14 at 18:19
  • 3
    $\begingroup$ My fear is that the OP has little sense about what he is asking. $\endgroup$ – wolfies Nov 29 '14 at 1:24
5
$\begingroup$

[Update: Corrected error in and improved first method.]

Given the upvotes on the question, there seems to be, or to have been, some interest in this question. As has been mentioned, based on a continuous probability model, the chances of getting a 3. is nil. But RandomReal is a discrete approximation so there is some small, positive probability of getting a particular number.

I can read the question in two ways, to avoid getting the particular number 3. or to avoid getting numbers Equal to the number 3. These mean different things in Mathematica.

No 3.'s

One way to avoid getting a 3. is to replace them like this:

pos3 = Compile[{{x, _Real, 1}}, Position[x, 3.], RuntimeOptions -> "Speed"];
l = ReplacePart[#, pos3[#] :> NestWhile[RandomReal[{3, 6}] &, 3., pos3[{#}] =!= {} &]] &@
      RandomReal[{-3, 6}, 1000]

It turns out a sure way to match exactly the number 3. is to use a compiled function with the RuntimeOptions setting "CompareWithTolerance" -> False. This is achieved with RuntimeOptions -> "Speed" along with some other speed-ups, like no overflow checking. See the section at the end for a little hitch with the replacement /. 3. -> .., which matches more than 3. and cannot be used to strictly eliminate only threes.

None equal to 3.

To avoid getting a number equal to 3., we need to understand how Equal works on approximate real numbers. It compares with a tolerance that is determined by Internal`$EqualTolerance, although testing shows that is not exactly true. It seems to depend on the bits in the binary representation of the floating-point number. The thresholds for equality with 3. should be about

3. ± 0.5 * (3 * 10^Internal`$EqualTolerance * $MachineEpsilon)

Where the 0.5 comes from the approximate boundary for rounding toward 3. We can solve for the exact boundaries, by first approximating them with FindRoot:

thresholdPlus = x /. FindRoot[
   Boole[3. == 3. + x*3*10^Internal`$EqualTolerance*$MachineEpsilon] - 1/2,
   {x, 0.5, 0.505}]
thresholdMinus = x /. FindRoot[
   Boole[3. == 3. - x*3*10^Internal`$EqualTolerance*$MachineEpsilon] - 1/2,
   {x, 0.49, 0.5}]

FindRoot::brmp: The root has been bracketed as closely as possible with machine precision...

(*  0.502604  *)

FindRoot::brmp: The root has been bracketed as closely as possible with machine precision...

(*  0.497396  *)

We can adjust these thresholds manually until we find the last point at which Equal returns True:

3. == 3. + (1 - 2 $MachineEpsilon) thresholdPlus *
    3 * 10^Internal`$EqualTolerance * $MachineEpsilon
3. == 3. - (1 - $MachineEpsilon) thresholdMinus * 
    3 * 10^Internal`$EqualTolerance * $MachineEpsilon
(*
  True
  True
*)

And reset the thresholds:

thresholdPlus = (1 - 2 $MachineEpsilon) thresholdPlus; 
thresholdMinus = (1 - $MachineEpsilon) thresholdMinus;

Now construct a function to use to pick out the 3.'s in a list:

equal3[x_] := 
  UnitStep[x - (3. - thresholdMinus*3*10^Internal`$EqualTolerance*$MachineEpsilon)] *
   UnitStep[(3. + thresholdPlus*3*10^Internal`$EqualTolerance*$MachineEpsilon) - x];

Then finally we can construct our list of random numbers like this:

l = FixedPoint[
 ReplacePart[#, SparseArray[equal3@#]["NonzeroPositions"] :> RandomReal[{-3, 6}]] &,
 RandomReal[{-3, 6}, 1000]]

It takes only about 4-6 times longer than a simple RandomReal, not counting the time it takes to figure out how to set up the code, which is significant.

Performance

TableForm[
 Table[
  {"10"^n,
   RandomReal[{-3, 6}, 10^n]; // RepeatedTiming // First,
   ReplacePart[#, 
         pos3[#] :> NestWhile[RandomReal[{3, 6}] &, 3., pos3[{#}] =!= {} &]
         ] &@RandomReal[{-3, 6}, 10^n]; // RepeatedTiming // First,
   FixedPoint[
       ReplacePart[#, 
         SparseArray[equal3@#]["NonzeroPositions"] :> 
          RandomReal[{-3, 6}]] &,
       RandomReal[{-3, 6}, 10^n]]; // RepeatedTiming // First,
   RandomVariate[
       MixtureDistribution[{6 - 2 $MachineEpsilon, 
         3 - 2 $MachineEpsilon}, {UniformDistribution[{-3, 
           3. - 2 $MachineEpsilon}], 
         UniformDistribution[{3. + 2 $MachineEpsilon, 6}]}], 10^n]; //
      RepeatedTiming // First},
  {n, 2, 7}],
 TableHeadings -> {None, {Length, RandomReal, ReplacePart, SparseArray, "Mixture"}}]

Mathematica graphics

The first method I gave (called ReplacePart in the table) is a bit faster than the SparseArray method, especially on shorter lists. Both add some noticeable overhead to a straight RandomReal call. In the OP's use-case (length 1000), ReplacePart is a bit under 3 times slower. Of course, as others have mentioned, it is highly unlikely that any replacement will ever be done. Note, also, how slow MixtureDistribution is.

What matches 3.? -- Caveat to the first solution

Interestingly, 3. matches two numbers (which I didn't know before):

{(1. - $MachineEpsilon) 3., 3. - 3*0.51 $MachineEpsilon, 3., 
  3. + 3*0.51 $MachineEpsilon, (1. + $MachineEpsilon) 3.} // FullForm
(*
  {2.999999999999999`,
   2.9999999999999996`,
   3.`,
   3.0000000000000004`, 
   3.000000000000001`}
*)

{(1. - $MachineEpsilon) 3., 3. - 3*0.51 $MachineEpsilon, 3., 
  3. + 3*0.51 $MachineEpsilon, (1. + $MachineEpsilon) 3.} /. 3. -> Infinity
(*  {3., ∞, ∞, 3., 3.}  *)

But if the OP is trying to avoid division by zero in, say, 1/(x-3), then the lower number would not need to be replaced:

1/((3. - 3*0.51 $MachineEpsilon) - 3)
(*  -2.2518*10^15  *)
$\endgroup$
3
$\begingroup$
data = DeleteCases[Range[-3, 6], 3];
l = RandomChoice[data, {1000}]
$\endgroup$
  • 4
    $\begingroup$ thanks for your answer but I said I want 1000 real number between -3 and 6 with 3 not in them.your answer outputs integers. $\endgroup$ – user2838619 Nov 27 '14 at 18:52
  • 2
    $\begingroup$ @user2838619 If you use RandomReal there are no integers anyway. So, what would you expect to obtain. $\endgroup$ – Alexei Boulbitch Nov 28 '14 at 9:20
1
$\begingroup$

One way to do this is to "redraw" a random number whenever a draw selects a number you wish to eliminate. Here is a compact and recursive way to do this :

random := With[{x = RandomReal[{-3, 6}]}, If[x == 3., random, x]];
Table[random, 1000]

(theoretically the code is not guaranteed to stop, although infinite running time has probability = 0).

To check that this works you can try

r := With[{x = RandomInteger[{-3, 6}]}, If[x == 3, r, x]];
data = Table[r, 1000];
Select[data, # == 3 &]

that returns

{}
$\endgroup$

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.