Segment 39

From Computational Statistics (CSE383M and CS395T)
Jump to navigation Jump to search

Calculation Problems

1. Calculate <math>\alpha(x_1, x_2)</math> for the given problem.

<math>\alpha(x_1, x_2)= min \left(1, \frac{\pi(x_{2c}) q(x_i|x_{2c})}{\pi(x_{1}) q(x_{2c}|x_{1})} \right)</math>

Acceptance Probability
x_2</math> 1 2 3 4 5
1 0 1 0 0 0
2 0.5 0 0.5 0 0
3 0 0.5 0 0.5 0
4 0 0 0.5 0 0.5
5 0 0 0 1 0

2.

To Think About

Notes: <math> \alpha = {O,P,G}</math> and <math> P_{\alpha} = \frac{W_\alpha N_\alpha}{\sum_\alpha W_\alpha N_\alpha} </math>


1. Gibb's Sampler:

def BallType(ball):
    if ball < 7: return 0
    if ball < 10: return 1
    if ball < 15: return 2
    raise Exception

#1-7 Orange 8-10 Purple 11-15 Green
a = [4,8,2,11,12,6]
def ProbGiven(draw):
    w = [6,4,3]
    urn = [7,3,5]
    probability = 1.0
    for i in draw:
        index = BallType(i)
        denom = sum(w[j]*urn[j] for j in range(len(w)))
        probability *= w[index]*urn[index]*1.0/denom
        urn[index] -=1
    return probability 

prob = ProbGiven(a)
for mmm in range(1):
    b = []
    for i in range(len(a)):
        replacements = set(range(len(a)))-set(a)
        replacements.add(a[i])
        replacements = list(replacements)
        print replacements
        p = []
        for j in replacements:
            a[i] = j
            p.append(ProbGiven(a))
        p = [j/sum(p) for j in p]
        print p
        choice = numpy.random.multinomial(1, p)
        for j,k in enumerate(choice):
            if k == 1:
                choice = j
        b.append(replacements[choice])
    a = b