PBRT 13.4 - computing the weight in Metropolis sampling

Practical and theoretical implementation discussion.
Iron
Posts: 1
Joined: Sat Jun 15, 2013 2:35 am

PBRT 13.4 - computing the weight in Metropolis sampling

Postby Iron » Mon Jun 17, 2013 11:02 am

I'm trying to implement the pseudocode on page 654 using equation 13.11 on page 657, but I don't understand how the 'weight' is computed.

Here is the pseudocode on page 653, and my implementation of it in sage, which, as far as I can tell, gives correct results:

Code: Select all

X = XO
for i = 1 to n
    X' = mutate(X)
    a = accept(X, X')
    if (random() < a)
       X = X'

Code: Select all

def F(x):
    """equation 13.11 on page 657"""
    if 0 <= x and x <= 1:
        return pow(x - 0.5, 2)
    else:
        return 0.00000001

samples = []
currentX = random()  #random starting point
n = 5000
for i in range(n):
    proposedX = mutate(currentX, 0, 1)
    a = min(1.0, F(currentX) / F(proposedX))
   
    if random() < a:
        currentX = proposedX
    samples.append( (currentX,  F(currentX)) )  #recording currentX with a value of F(currentX)


list_plot(samples)



The pseudocode on page 654, and my implementation in sage, which gives incorrect results:

Code: Select all

X = XO
for i = 1 to n
    X' = mutate(X)
    a = accept(X, X')
    record(X, (1-a) * weight)
    record(X', a * weight)
    if (random() < a)
       X = X'

Code: Select all

# ... as above
for i in range(n):
    proposedX = mutate(currentX, 0, 1)
    a = min(1.0, F(currentX) / F(proposedX))
   
    samples.append( (currentX, (1.0 - a) * F(currentX)) )
    samples.append( (proposedX, a * F(proposedX)) )
    if random() < a:
        currentX = proposedX

list_plot(samples)


So, how is the weight computed ?

Return to “General Development”

Who is online

Users browsing this forum: No registered users and 1 guest