3

Is there any build in function or shorter version of this? Assume I have 3 reactions time and each reaction occurs randomly. Generate rand=RandomReal[] if rand<0.2 then result is 1, if rand<0.5 then result is 2 and if rand<1 then result is 3. Thanks..

 reaction = {0.2, 0.5, 1};
 First@Flatten@Position[reaction, SelectFirst[reaction, RandomReal[] < # &]]
gwr
  • 13,452
  • 2
  • 47
  • 78
OkkesDulgerci
  • 10,716
  • 1
  • 19
  • 38

4 Answers4

7

Assuming you meant to only use one RandomReal call (per @Coolwater's comment), you can use RandomChoice:

RandomChoice[{.2, .3 ,.5} -> {1, 2, 3}]
RandomChoice[{.2, .3 ,.5} -> {1, 2, 3}, 10]

3

{3, 1, 3, 2, 3, 2, 3, 2, 3, 1}

Addendum

To mimic the OP code with RandomChoice, by changing the random number for each comparison, you just need to adjust the probabilities. 20% of the time, the first random number is less than .2. Of the other 80%, 50% of the time the second random number is less than .5, and 50% of the time it is greater than .5. So, the following should produce the same distribution as the OP function:

RandomChoice[{.2, .4, .4} -> {1, 2, 3}]

We can check by running the OP code multiple times:

reaction = {0.2, 0.5, 1};
data = Table[First@Flatten@Position[reaction, SelectFirst[reaction, RandomReal[] < # &]], {10^5}];

Counts[data]/10.^5

<|1 -> 0.20025, 2 -> 0.4018, 3 -> 0.39795|>

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
  • What I want is generate random number rand=RandomReal[] and look at its location in reaction and extract index of reaction. – OkkesDulgerci Dec 27 '17 at 18:45
  • @OkkesDulgerci I think that's exactly what my function does. 20% of the time, the random number is less than .2, 30% of the time the random number is between .2 and .5 and 50% of the time it is more than .5. – Carl Woll Dec 27 '17 at 18:53
  • RandomChoice[{.2, .3 ,.5} -> {1, 2, 3}] might be equivalent what I am doing. I will think about it. – OkkesDulgerci Dec 27 '17 at 18:54
  • When I replaced your code with mine it does not give desired result. – OkkesDulgerci Dec 27 '17 at 19:06
  • This looks right to me. Clarify why the result is not as desired. – george2079 Dec 27 '17 at 19:08
  • what if reaction is reaction={0.999998, 0.999999, 1.}; – OkkesDulgerci Dec 27 '17 at 19:09
  • 1
    @OkkesDulgerci - your code (generating a random number for each test done by SelectFirst) is equivalent to RandomChoice[{0.2, 0.4, 0.4} -> Range[3]] since the first test extracts ~20%, if the second test is needed it extracts ~50% of the remaining ~80% (i.e., ~40%), and the last test will extract all of the remaining. Using a single random number gives the result that Carl provides RandomChoice[{0.2, 0.3, 0.5} -> Range[3]] – Bob Hanlon Dec 27 '17 at 19:12
  • One can find complete code here https://mathematica.stackexchange.com/questions/158300/michaelis-menten-kinetics-using-gillespie-algorithm/158543#158543 – OkkesDulgerci Dec 27 '17 at 19:15
  • your question implies {.2,.5,1} is fixed. The generalization is RandomChoice[Differences@Prepend[reaction, 0] -> Range[3], 100] (assume the last element is always 1 ) – george2079 Dec 27 '17 at 19:17
  • I liked @CarlWoll solution but for some reason it does not work. – OkkesDulgerci Dec 27 '17 at 19:17
  • Yes last element in reaction is always 1 sorry for being unclear. – OkkesDulgerci Dec 27 '17 at 19:18
3
wd = WeightedData[Range[3], Differences[{0, ## & @@ reaction}]]; 

SeedRandom[1]
RandomChoice[wd]

3

RandomChoice[wd, 5]

{1,3,1,2,1}

Update: a shorter alternative using FirstPosition

SeedRandom[123]
FirstPosition[reaction - RandomReal[], _?Positive][[1]]

2

Note: This function uses a single random number while the alternative

FirstPosition[reaction, _?(#>RandomReal[]&)][[1]]

creates a new random number for every comparison until the condition is satisfied (as does the code in OP), and hence, it can give a different result 

SeedRandom[123]
FirstPosition[reaction, _?(#>RandomReal[]&)][[1]]

3

kglr
  • 394,356
  • 18
  • 477
  • 896
  • I'll go with this since it is more readable for me. First@FirstPosition[reaction,_?(RandomReal[]<#&)] – OkkesDulgerci Dec 27 '17 at 21:22
  • First@FirstPosition[reaction,_?(RandomReal[]<#&)] does not work in Do loop but this works perfectly FirstPosition[reaction - RandomReal[], _?Positive][[1]] see https://mathematica.stackexchange.com/questions/165073/birth-death-process-with-delay/172031#172031 – OkkesDulgerci Apr 28 '18 at 07:46
1

Another possibility is to make this a distribution and then use all the machinery that goes with that:

SeedRandom["December 27, 2017"];
distReaction = With[
    {
       reaction = {0.2, 0.5, 1}
    },
    EmpiricalDistribution[ Differences[{0} ~ Join ~ reaction ] -> Range[3] ]
]


RandomVariate[ distReaction, 5 ]

{2, 1, 3, 3, 3}

We could then ask for the probability of this specific result assuming that the results are i.i.d.:

Probability[ 
    {x1, x2, x3, x4, x5} == {2, 1, 3, 3, 3}, 
    Thread [ {x1, x2, x3, x4, x5} \[Distributed] distReaction ]
]

0.0075

gwr
  • 13,452
  • 2
  • 47
  • 78
1

The first response that came to mind, after reading the question, was 

BlockRandom[
 With[{rand = RandomReal[]},

  Which[
   rand <= 0.2, reaction[[1]],
   0.2 < rand <= 0.5, reaction[[2]],
   0.5 < rand <= 1, reaction[[3]]]], RandomSeeding -> 123456798]

where, obviously, reaction = {0.2, 0.5, 1}. By replacing reaction[[i]] with i it is easy to obtain the index of the reaction (not the time it takes).

I think that, this is as literal an implementation as it gets.

A couple of other approaches are

BlockRandom[
 With[{rand = RandomReal[]},

  Extract[
   reaction,
   First[Position[Thread[rand <= reaction], True]]]
  ], RandomSeeding -> 123456798]

and 

BlockRandom[
 With[{rand = RandomReal[]},

  Last[Flatten[Reap[
     Scan[
      If[rand <= #, Sow[#]; Return[]] &, reaction]]]]
  ], RandomSeeding -> 123456798]

In the former of the last two approaches, the index of the reaction instead of it's time can be obtained by replacing Extract[<>] with First[Flatten[Position[Thread[rand <= reaction], True]]].

The later code segment cannot be readily adjusted to obtain the index of the reaction instead of the time (at least not without doing something like First@Position[reaction, x] where x stands for the output of the last approach)

user42582
  • 4,195
  • 1
  • 10
  • 31