Massive Math Stuff: Probabilities

Showing posts with label Probabilities. Show all posts

Saturday, November 21, 2015

Coupon collector's problem

https://en.wikipedia.org/wiki/Coupon_collector%27s_problem
the coupon collector's problem describes the "collect all coupons and win" contests. It asks the following question: Suppose that there is an urn of n different coupons, from which coupons are being collected, equally likely, with replacement. What is the probability that more than tsample trials are needed to collect all n coupons? An alternative statement is: Given n coupons, how many coupons do you expect you need to draw with replacement before having drawn each coupon at least once? The mathematical analysis of the problem reveals that the expected number of trials needed grows as

\Theta(n\log(n))

.^[1] For example, when n = 50 it takes about 225^[2] trials to collect all 50 coupons.

赠券收集问题

赠券收集问题的特征是开始收集时，可以在短时间内收集多种不同的赠券，但最后数种则要花很长时间才能集齐。例如有50种赠券，在集齐49种以后要约多50次收集才能找到最后一张，所以赠券收集问题的答案t的期望值要比50要大得多。

http://www.cnblogs.com/lailailai/p/4886389.html

以一个例子来说明：有种食品每袋都会附赠一张卡片，共有10种卡片，每种卡片出现的概率相等。

首先来求当拥有i种卡片时，再拆开一袋得到不同卡片的概率Pi：

a. 当我们没有任何卡片时，随意拆开就能得到一种不同的卡片，即第一种得到的卡片P0 = 1

b. 当有了一种继续拆，此时得到新的种类的概率为P1 = 9 / 10
...当有了i种，此时得到新种类的概率为Pi = (10 - i) / 10

我们知道了Pi的值，那么对于拆了i包后，拆开发现新卡片的次数期望为E = Pi * Ni，由于只需要然其发生一次即可令E=1反过来求对应的Ni，Ni = 1 / Pi，即拆1 / Pi包时可以发现一种新卡片（相对于已经拥有i种卡片）。

那么发现全部卡片总共需要的包数 = N0 + N1 + N2 ... N9 = 10 / 10 + 10 / 9 + 10 / 8 ... 10 / 1 = 29.289682539682538

即按期望来，买30包可以集齐卡片。

下面给出了卡片等概率出现时，集齐卡片的期望次数：

1 1.0
2 3.0
3 5.5
4 8.33333333333
5 11.4166666667
6 14.7
7 18.15
8 21.7428571429
9 25.4607142857
10 29.2896825397
11 33.2186507937
12 37.2385281385
13 41.3417388167
14 45.5218725719
15 49.7734348984
16 54.0916638917
17 58.4723928849
18 62.9119454075
19 67.4070534857
20 71.9547931429
21 76.5525328
22 81.1978915048
23 85.888704755
24 90.6229962661
25 95.3989544438
26 100.214912622
27 105.069332338
28 109.960789091
29 114.88796013
30 119.849613928
31 124.844601059
32 129.871846254
33 134.930341449
34 140.019139675
35 145.137349666
36 150.284131085
37 155.458690281
38 160.660276505
39 165.888178519
40 171.141721557
41 176.420264596
42 181.723197879
43 187.049940686
44 192.399939306
45 197.7726652
46 203.167613315
47 208.584300561
48 214.022264403
49 219.481061578
50 224.960266916
51 230.459472255
52 235.978285436
53 241.516329387
54 247.073241262
55 252.648671656
56 258.242283868
57 263.853753223
58 269.482766437
59 275.129021031
60 280.792224777
61 286.47209519
62 292.168359046
63 297.880751933
64 303.609017837
65 309.352908741
66 315.11218426
67 320.886611294
68 326.675963702
69 332.480021991
70 338.298573035
71 344.131409792
72 349.978331057
73 355.839141211
74 361.713649994
75 367.601672291
76 373.503027922
77 379.417541447
78 385.345041986
79 391.285363037
80 397.238342316
81 403.203821595
82 409.181646553
83 415.171666632
84 421.173734905
85 427.18770794
86 433.21344568
87 439.250811328
88 445.299671228
89 451.359894765
90 457.431354256
91 463.513924859
92 469.607484473
93 475.711913652
94 481.827095519
95 487.952915684
96 494.089262165
97 500.236025313
98 506.393097739
99 512.560374246
100 518.737751764
101 524.925129282
102 531.122407789
103 537.329490219
104 543.546281386
105 549.772687938
106 556.008618299
107 562.253982622
108 568.50869274
109 574.772662118
110 581.045805807
111 587.328040405
112 593.619284012
113 599.919456191
114 606.228477927
115 612.546271593
116 618.872760911
117 625.207870919
118 631.551527936
119 637.903659528

E = n / n + n / (n-1) + n/ (n-2) ... n / 1

当卡片出现的概率不是相等的时候，Pi的计算过程就麻烦许多

Tuesday, November 10, 2015

Puzzles, Maths and Algorithms: Unbiased Coin Tossing

Puzzles, Maths and Algorithms: Unbiased Coin Tossing
Given a biased coin, with probability of Heads equal to x. How to do unbiased coin tossing?

Lets define an event, E as tossing the biased coin twice. The possible outcomes with probabilities is as follows

P(h,h) = x^2
P(h,t) = x(1-x)
P(t,h) = x(1-x)
P(t,t) = (1-x)^2

The event h,t or t,h are equi-likely, without any bias we can call that if Event h,t occurs it means head, t,h means tails but if h,h or t,t occurs we repeat the experiment.

Try to find the expected number of coin toss that would be required to call heads or tails?

Let expected coin toss be e.

Probability that we get outcome in 1st event is 2x(1-x)
Total number of coin toss would be 2
Probability that we get no outcome in 1st event is [1 - 2x(1-x)]
Total number of coin toss would be 2 + e (we wasted 2 coin toss and still we expect e)

==> e = 2 * 2x(1-x) + (2+e) * [1 - 2x(1-x)]
==> e = 4x(1-x) + 2 - 4x(1-x) + e - 2ex(1-x)
==> e = 1/ [x(1-x)]

so for x = 1/2, e = 4
for x=2/3, e = 4.5
for x=1, e = INF (as expected, because all we get is chain of h h h h h)

What is the probability that there wont be any outcome in e coin tosses (expected outcomes)?

p = prob of e heads + prob of e tail
==> p = x^e + (1-x)^e

To find a bound on this p in polynomial terms can be done by using binomial expansion and Newtonian series is out of the scope of this blog. But some number crunching is.

For x = 1/2, e = 4, p = 0.125
For x = 2/3, e = 4.5, p = 0.168
For x = 3/4, e = 5.33, p = 0.216
For x = 0.99, e = 101, p = 0.362

This tells us that probability that outcome comes is quite good even for high coin biases.

How can we improve on the expected number of coin tosses?

In above method, we say that HT or TH terminates experiment. and continue the experiment a fresh when the outcome is HH or TT.

We can further combine outcome of two such events to increase the probability of outcome e.g. say HH TT => heads and TT HH means tails.

How much expected coin tosses are we doing in this case?

Read full article from Puzzles, Maths and Algorithms: Unbiased Coin Tossing

Puzzles, Maths and Algorithms: Stick Breaking

Puzzles, Maths and Algorithms: Stick Breaking
Problem Set 1: Given a stick of length 1 meter. It is cut into two parts randomly.

Problem 1a: What is the expected length of the smaller stick?
Problem 1b: What is the expected length of the bigger stick?

Problem Set 2: Given a stick of length 1 meter. It is cut into three parts randomly.

Problem 2a: What is the expected length of the smallest stick?
Problem 2b: What is the expected length of the largest stick?
Problem 2c: What is the expected length of the stick with second smallest (or second largest) length?
Problem 2d: What is the expected length of the middle stick (stick that is between the two cuts)?
Problem 2e: What is the probability that the three sticks would form a triangle?
Problem 2f: What is the probability that the three sticks would form a right angle triangle?
Problem 2g: What is the probability that the three sticks would form an equilateral triangle?

Problem Set 3: Given a stick of length 1 meter. It is cut into two parts randomly. Then bigger stick is picked and is again cut into two parts randomly.

Problem 3a: What is the expected length of the smallest stick?
Problem 3b: What is the expected length of the largest stick?
Problem 3c: What is the expected length of the stick with second smallest (or second largest) length?
Problem 3d: What is the expected length of the middle stick (stick that is between the two cuts)?
Problem 3e: What is the probability that the three sticks would form a triangle?
Problem 3f: What is the probability that the three sticks would form a right angle triangle?
Problem 3g: What is the probability that the three sticks would form an equilateral triangle?

Generalized Cut Problem: Find the solution to problem 2a-2c, 3a-3c for n-1 cuts.

Let us work out solution for 1a. Smaller stick has length x which is uniformly distributed in the range [0,0.5]. The probability distribution of x is p(x) = 1 / [1/2-0] = 2. Yes it is a pdf, ∫ ₀^0.5p(x)dx = 1. The expected length of the smaller stick is E[x] = ∫ x p(x) dx = ∫ ₀^0.52xdx = 1/4.

2a) Let us work out solution for 2a. Let two cuts are at distance x,y from one end. Let x < y. So the three sticks are of length x, y-x, 1-y. Let z = min(x,y-x,1-y) represent the minimum length of the stick. Let Pr(z>a) indicate the probability that min length stick has length greater than a.

=> Pr(z>a) = 2 * Pr(x>a, y-x>a, 1-y>a) = Pr(x>a, a+x < y < 1-a)

Note multiplication by 2 is for the other case (x > y). Since a+x< y < 1-a, we get x < 1-2a. Also setting x =a, we get a < 1/3 (which is obvious).

=> Pr(z > a) = 2 ∫ _a^1-2a ∫ _a+x^1-adxdy = (1-3a)².

So we get Pr(z<=a) = 1-Pr(z>a). And Pr(z=a) = derivative of Pr(z<=a) at a = 6(1-3a).

=> E[a] = ∫ ₀^1/36a(1-3a)da = 1/9

2b) We can do a similar derivation for z = max(x,y-x,1-y) to get the expected max length. It would be slightly more complicated. Unlike above, we compute here Pr(z< a)

=> Pr(z < a) = Pr(x < a, y-x < a, 1-y < a) = Pr( 1-2a < x < a, 1-a < y < x + a).

We need to be careful here are respect constraints like x > 0, x < y, y < 1. To ensure these constraints are respected. Pr(z < a) is broken into three cases, Case 1: a < 1/2. Case 2: a >= 1/2, 0 < x < 1-a, Case 3: a >= 1/2, 1-a < x < a. Then we have to take differential of obtained Pr(z < a) to get Pr(z=a). [Dont mix terms from the three cases].

Integrate the first case from [1/3 to 1/2]. Second and third case from [1/2 to 1]. We will see the desired output as 1/2+1/9=11/18. (Don't forget to multiply by 2 to cover for x > y case).

2c) It would be easy to obtain: 1 - (2a) - (2b) = 1 - 1/9 - 1/2 - 1/9 = 1/2 - 2/9 = 5/18.

2d) Since there are no constraints, expected length of middle cut is 1/3

2e) For triangle, we require satisfaction of three inequalities
x + y - x > 1-y,
x + 1- y > y - x,
y - x + 1 - y > x
Apart from the implicit x < y inequality.

=> 1/2 < y < x + 1/2, x < 1/2

So x follows a pdf of P(x) = 2. y satisfied the range with probability x. So the probability of triangle formation is ∫ ₀^1/22 x dx = 1/4.

2f) There are discrete stick lengths, e.g. (3/12, 4/12, 5/12) that satisfy this property. There are countably infinite, yet the area encapsulated by them is 0.

Problem set 3 is slightly more complicated. Let us work out 3a).

Let the first cut is at length x from shorter side (x < 1/2). Now from the stick of length (1-x), we make a cut at random of length y. Let y < (1-x)/2. If x is the shortest stick then x < y => x < 1/3.

Case 1: 1/3 < x < 1/2, 0 < y < (1-x)/2 => y is the shortest stick

Case 2: 0 < x < 1/3
Case 2a: 0 < y < x => y is the shortest stick
Case 2b: x < y < (1-x)/2 => x is the shortest stick

We also note that P(x) = 2 (a uniform distribution in [0,1/2]), P(y) = 2/(1-x) (a uniform distribution in [0, (1-x)/2]).

So the Expected min length = ∫ _1/3^1/2 ∫ ₀^(1-x)/2 4y/(1-x) dx dy + ∫ ₀^1/3 ∫ ₀^x 4y/(1-x) dx dy +
∫ ₀^1/3∫ _x^(1-x)/24x/(1-x) dx dy = 15/16 - 2 log(3/2)

For Generalized problem read David and Nagaraja's Order Statistics (pp. 133-135, and p. 153)

The solutions through Monte-Carlo simulations. Following python code solves problem set 2 and 3.

from random import random
def prob2():
    niter = 100000
    n, ntri, nrtri = 0.0, 0, 0
    lsmall, llarge, lmid, lcenter = 0, 0, 0, 0

    while n < niter:
        n += 1
        p = random()
        q = random()
        x = min(p, q)
        y = max(p, q)

        l1 = x
        l2 = y-x
        l3 = 1-y

        mn = min(l1, l2, l3)
        mx = max(l1, l2, l3)
        md = 1 - mn - mx
        lsmall += mn
        llarge += mx
        lmid += md
        lcenter += y

        if l1 + l2 > l3 and l1 + l3 > l2 and l2 + l3 > l1:
            ntri += 1

        if mn**2 + md**2 == mx**2:
            nrtri += 1

    print 'min len:', lsmall/n
    print 'max len:', llarge/n
    print 'mid len:', lmid/n
    print 'center stick len', lcenter/n
    print 'triangle prob:', ntri/float(n)
    print 'right triangle prob:', nrtri/float(n)

def prob3():
    niter = 100000
    n, ntri, nrtri = 0.0, 0, 0
    lsmall, llarge, lmid, lcenter = 0, 0, 0, 0

    while n < niter:
        n += 1
        p = random()
        q = random()

        l1 = 0.5 * (1-p)
        l2 = 0.25 * (1+p) * (1-q)
        l3 = 0.25 * (1+p) * (1+q)

        mn = min(l1, l2, l3)
        mx = max(l1, l2, l3)
        md = 1 - mn - mx
        lsmall += mn
        llarge += mx
        lmid += md
        lcenter += (l2+l3)/2

        if l1 + l2 > l3 and l1 + l3 > l2 and l2 + l3 > l1:
            ntri += 1

        if mn**2 + md**2 == mx**2:
            nrtri += 1

    print 'min len:', lsmall/n
    print 'max len:', llarge/n
    print 'mid len:', lmid/n
    print 'center stick len', lcenter/n
    print 'triangle prob:', ntri/float(n)
    print 'right triangle prob:', nrtri/float(n)

Read full article from Puzzles, Maths and Algorithms: Stick Breaking