Pólya urn flip and roll

13

4

Problem statement

Pólya is playing about with his urn again and he wants you to help him calculate some probabilities.

In this urn experiment Pólya has an urn which initially contains 1 red and 1 blue bead.

For every iteration, he reaches in and retrieves a bead, then inspects the colour and places the bead back in the urn.

He then flips a fair coin, if the coin lands heads he will insert a fair 6 sided die roll amount of the same coloured bead into the urn, if it lands tails he will remove half the number of the same colored bead from the urn (Using integer division - so if the number of beads of the selected colour is odd he will remove (c-1)/2 where c is the number of beads of that colour)

Given an integer n ≥ 0 and a decimal r > 0, give the probability to 2 decimal places that the ratio between the colours of beads after n iterations is greater than or equal to r in the shortest number of bytes.

An example set of iterations:

Let (x, y) define the urn such that it contains x red beads and y blue beads.

Iteration    Urn       Ratio
0            (1,1)     1
1            (5,1)     5        //Red bead retrieved, coin flip heads, die roll 4
2            (5,1)     5        //Blue bead retrieved, coin flip tails
3            (3,1)     3        //Red bead retrieved, coin flip tails
4            (3,4)     1.333... //Blue bead retrieved, coin flip heads, die roll 3

As can be seen the Ratio r is always ≥ 1 (so it's the greater of red or blue divided by the lesser)

Test cases:

Let F(n, r) define application of the function for n iterations and a ratio of r

F(0,5) = 0.00
F(1,2) = 0.50
F(1,3) = 0.42
F(5,5) = 0.28
F(10,4) = 0.31
F(40,6.25) = 0.14

This is code golf, so the shortest solution in bytes wins.

Expired Data

Posted 2019-03-15T15:15:33.260

Reputation: 3 129

I feel like there is a formula for this... – Embodiment of Ignorance – 2019-03-15T18:10:42.853

Something to do with beta binomials maybe, but it might be longer to write that out – Expired Data – 2019-03-15T18:13:54.027

depends on the language; R and Mathematica might be able to do it efficiently. – Giuseppe – 2019-03-15T20:45:12.777

Answers

6

JavaScript (ES7),  145 ... 129 124  123 bytes

Takes input as (r)(n). This is a naive solution that actually performs the entire simulation.

r=>g=(n,B=s=0,R=0,h=d=>++d<7?h(d,[0,d].map(b=>g(n,B/-~!!b,R/-~!b)&g(n,B+b,R+d-b))):s/24**-~n)=>n--?h``:s+=~B<=r*~R|~R<=r*~B

Try it online!

Too slow for the last 2 test cases.

Commented

r =>                    // r = target ratio
g = (                   // g is a recursive function taking:
  n,                    //   n = number of iterations
  B =                   //   B = number of blue beads, minus 1
  s = 0,                //   s = number of times the target ratio was reached
  R = 0,                //   R = number of red beads, minus 1
  h = d =>              //   h = recursive function taking d = 6-sided die value
    ++d < 7 ?           // increment d; if d is less than or equal to 6:
      h(                //   do a recursive call to h:
        d,              //     using the new value of d
        [0, d].map(b => //     for b = 0 and b = d:
          g(            //       do a first recursive call to g:
            n,          //         leave n unchanged
            B / -~!!b,  //         divide B by 2 if b is not equal to 0
            R / -~!b    //         divide R by 2 if b is equal to 0
          ) & g(        //       do a second recursive call to g:
            n,          //         leave n unchanged
            B + b,      //         add b blue beads
            R + d - b   //         add d - b red beads
          )             //       end of recursive calls to g
        )               //     end of map()
      )                 //   end of recursive call to h
    :                   // else (d > 6):
      s / 24 ** -~n     //   stop recursion and return s / (24 ** (n + 1))
) =>                    // body of g:
  n-- ?                 //   decrement n; if n was not equal to 0:
    h``                 //     invoke h with d = [''] (coerced to 0)
  :                     //   else:
    s +=                //     increment s if:
      ~B <= r * ~R |    //       either (-B-1) <= r*(-R-1), i.e. (B+1)/(R+1) >= r
      ~R <= r * ~B      //       or     (-R-1) <= r*(-B-1), i.e. (R+1)/(B+1) >= r

Arnauld

Posted 2019-03-15T15:15:33.260

Reputation: 111 334

I really like this answer, I found that in order to solve the later test cases I needed to add code to merge the same ratio probabilities. So I'm not surprised it's too slow – Expired Data – 2019-03-15T20:54:54.430

0

Wolfram Language (Mathematica), 133 bytes

Length@Select[#2/#&@@@(r=Nest[Sort/@Flatten[Table[{{k,0}+#,⌈#/{1,2}⌉},{k,6}]&/@Join[#,Reverse/@#],2]&,{{1,1}},x]),#>=y&]/Length@r


By using N[F,2],the output should print 2 decimal digits but TIO prints more...

Try it online!

J42161217

Posted 2019-03-15T15:15:33.260

Reputation: 15 931