From Computational Statistics (CSE383M and CS395T)
< User:Jzhang
Revision as of 21:37, 18 January 2013 by Jzhang (talk | contribs) (Activity 2)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Solutions for segment 2.

Skilled problem

Problem 1: If the knight had captured a Gnome instead of a Troll, what would his chances be of crossing safely?


<math>H_1</math> be the hypothesis that he is on bridge without troll
<math>H_2</math> be that he is on bridge with one troll
<math>H_3</math> be that he is on bridge with two trolls

The probability of crossing safely(or on the bridge without troll) is:

<math> P(H_1|G) = \frac{P(G|H_1)*P(H_1)}{P(G)} = \frac{P(G|H_1)*P(H_1)}{P(GH_1)+P(GH_2)+P(GH_3)} = \frac{1*\frac{3}{5}}{1*\frac{3}{5}+\frac{4}{5}*\frac{1}{5}+\frac{3}{5}*\frac{1}{5}} = \frac{15}{22}</math>

Problem 2: Suppose that we have two identical boxes, A and B. A contains 5 red balls and 3 blue balls. B contains 2 red balls and 4 blue balls. A box is selected at random and exactly one ball is drawn from the box. What is the probability that it is blue? If it is blue, what is the probability that it came from box B?

The probability of be blue ball is: <math> P(b) = P(bA)+P(bB) = \frac{1}{2}*\frac{3}{8}+\frac{1}{2}*\frac{4}{6} = \frac{25}{48}</math>

The probability of being box B given blue ball is: <math> P(B|b) = \frac{P(b|B)*P(B)}{P(b)} = \frac{16}{25}</math>

Thought problem

All the work is done in team with Sean Trettel

Activity 1

  • Simulate the Knight/Troll/Gnome problem 100,000 times.
  • Plot (fraction of safe crossings so far) vs. (number of simulated trials so far) to confirm that this fraction converges to the probability calculated in the segment.

The idea is simulating the whole world with a list of bridges, each bridge is represented with a list, there're five species on each bridge, which are represented with 5 numbers in the list, 1 means Gnome, -1 means Troll.

In each run of the simulation, a random integer number between 0~4 will determine which bridge the knight is on, and another random integer number between 0~4 will determine which specie the knight capture. If the knight capture a Troll, the number of capturing Troll is incremented, and whether he can cross safely is calculated (If the rest four species are Gnomes, it's safe) by summing rest of the list, if it's equal to 4, there're 4 Gnomes left, it's safe, and number that he can cross safely is incremented.

Here's the python code:

import pylab
import numpy
import random
b1 = [1,1,1,1,1]
b2 = [-1,1,1,1,1]
b3 = [-1,-1,1,1,1]

world = [b1,b1,b1,b2,b3]

cTroll = 0
cSafe = 0
frac_list = []

for i in range (10**5):
    # choose the world
    bc = world[random.randint(0,4)]
    # capture one random specie
    capture = bc[random.randint(0,4)]
    if capture < 0: # capture the troll
        if sum(bc)-capture == 4: # cross safely
    # record fraction at every step
    if not cTroll:

# plot the fraction verse the simulation time
x = numpy.arange(0,10**5,1)

pylab.xlabel('simulation times')

print cSafe
print cTroll
print cSafe*1.0/cTroll

The calculated fraction is 0.329046526868.
The plot is shown here:

Frac vs simulation.png

Activity 2

  • You are an oracle that, when asked, says "yes“ with probability ½ and "no" with probability ½. How do you do this using only a coin that comes up heads with unknown but constant probability P?
  • Simulate your scheme to confirm that it works.

The initial idea is roughly estimate the p first, and deduce some hypothesis that has an probability of 1/2.

The p can be estimated by flipping the coins a couple of thousand times and calculate the fraction of the times with head facing up. And then, we can create an artifical case, such as through coin n times, the chance that the first n-1 time the coin is tail and the nth time it's head is following the geometrical distribution:

<math> P = (1-p)^np</math>
Since we has the estimation of p, we can solve the n in the equation <math> P = 1/2 </math>

Now we know if we flip the coin n times, the chance of the first n-1 time being tail and last time being head is 1/2, so if this happened, we say yes, otherwise we say know.

The problem about this idea is that the n we solve is not necessary to be a number larger than 1 or even positive, for example if p is 0.9, the n we solve is -14.27, there's no way we can flip -14 times.

This fits sorts to the idea of transformation.

The implementation of the really elegant way of solving this problem with rejection is here:

import random

# given probability

# flip coins twice each round, only count the times of tail/head or head/tail


for i in range(10**6):
    # flip coin for the first time
    if random.random() < p:
        r1 = 1
        r1 = 0
    # flip coin for the second time
    if random.random() < p:
        r2 = 1
        r2 = 0
    # head/tail or tail/head
    if r1 + r2 == 1:
        if r1:

print p
print yes
print no
print yes*1.0/(yes+no)