## 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

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 = XOfor 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.00000001samples = []currentX = random()  #random starting pointn = 5000for 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 = XOfor 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 abovefor 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 = proposedXlist_plot(samples)`

So, how is the weight computed ?