I am going through Bayesian Melding paper by Poole and Raftery (2000). One of the ideas of the paper is demonstrated by Example 3.5, where there are three uniform distributions considered for $$X\sim U[2,4], Y\sim U[6,9], Z\sim U[0,5]$$. These distributions can be thought as a prior distribution of these variables. However, there is one constraint on the variables: $$Z= X/Y~.$$ Given this constraint and the prior distributions, the idea is obtain posterior distribution for $$Z, X, Y$$ through a method called pooling. In this method the following would be done to implement Example 3.5 numerically:

1. Sample say $$n=10^5$$ samples from $$X,Y$$:
X = np.random.uniform(2, 4, size=100000)
Y = np.random.uniform(6, 9, size=100000)

1. Calculate $$Z$$ from the generated sample of $$X, Y$$:
Z= X/Y

1. Calculate the probability density of the generated $$Z$$ through non-parametric method. In this case we can use Gaussian kernel to calculate pdf of generated $$Z$$, $$\hat{p}_Z$$.
def density_estimator(Z, h=0.1):
'''
phi: array of Z values
return: gaussian kernel smoothed point estimate of density of each Z
'''
density_return = np.copy(Z)
for i in range(len(Z)):
Z_curr = Z[i]
density_return[i] = np.sum(np.exp(-(Z_curr-Z)**2/(2*h**2)))
return 1./len(Z)*1./(2*np.pi*h**2)**0.5*density_return

estimated_density = density_estimator(Z)

1. Calculate importance sampling weights $$w_k=\sqrt{\frac{p_Z(Z)}{\hat{p}_Z(Z)}}$$. Because the true distribution of $$Z$$ is uniform, and the span of generated $$Z$$ is within the range covered by $$Z$$, the weights are simply proportional to:
weights = 1./estimated_density**0.5
#normalise
weights /= np.sum(weights)

1. Sample $$l$$ values from the discrete distributions (obtained in steps 1 and 2) with values $$X, Y$$ and $$Z$$ and probabilities proportional to weights. This is the step I am confused about. I did the following:
n = 10000
smp = np.copy(Z[:n]) #initialise the size of posterior sample of Z

def importance_sampling(smp, weights, Z):
j, k = 0, 0
while j < len(smp):
i = randint(0, weights.size)  # Random position of w
if random() <= weights[i]:
smp[j] = Z[i]
j += 1
k += 1
return smp
sample_Z = importance_sampling(smp, weights, Z)


If I do the following, the posterior distribution of $$Z$$ I obtain looks very different from Figure 3c of the paper. Any idea on where I am missing things?

• $X\sim U[2,4], Y\sim U[6,9], Z\sim U[0,5],Z= X/Y$ do not seem consistent Jun 6, 2021 at 13:26
• @Henry That's the whole premise of this problem, if you have prior distributions which follow certain relation, the idea is to find a compromise between the prior knowledge and the relation to obtain the posterior distributions. You can read the whole idea in the paper link. Jun 6, 2021 at 20:54