0
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ $X\sim U[2,4], Y\sim U[6,9], Z\sim U[0,5],Z= X/Y$ do not seem consistent $\endgroup$
    – Henry
    Jun 6, 2021 at 13:26
  • $\begingroup$ @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. $\endgroup$
    – titanium
    Jun 6, 2021 at 20:54

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.