3
$\begingroup$

At the moment I'm working with an algorithm that calculates probabilities $Z$ for a finite sequence of data points $X= \left\{ x_1,x_2,...,x_n \right\}$ using an expression like:
$p(Z_i| X_i,\theta) = \frac{1}{\phi}f(X_i,\theta)$,
where $\phi$ is a normalising constant such that $\sum_n p(Z_n|X_n) = 1$. (i.e. $f(X_i,\theta)$ gives us a number, we have a sequence of numbers and we normalise over all of them by summing and dividing). The next part of the algorithm chooses one of the $x$'s with weighted probability.
At the moment, it is possible to enumerate the sample space and calculate $\phi$, but is it possible to have some kind of MCMC approach, to avoid having to do this? So sampling a few positions, and evaluating $p(Z_i|X_i)$ only at those positions? I can see how MCMC works in the continuous case, but not sure if it's applicable here.
If it makes any difference, the values of $f(X_i,\theta)$ can have a lot of variation (e.g. 0.08, 0.70, 0.26, 64.00, 0.79, 7.11, 0.01, etc.) so areas of high probability can be separated by regions of low probability.

Thanks in advance.

$\endgroup$
  • $\begingroup$ I don't see how this problem is solvable without making additional assumptions about $f$. For instance, suppose the sample space has $10^{100}$ elements, the values of $f$ just happen to sum to unity, but exactly one value of $f$ equals $1 - 10^{-100}$. In other words, most of the probability is concentrated at an extremely small subset. Even a long MCMC run has little chance of encountering such a subset and so its results will be grievously wrong. $\endgroup$ – whuber Dec 5 '12 at 16:00
  • $\begingroup$ Is your end objective to generate a random sample of the $x_i$ from the distribution $f$ without having to calculate the normalizing constant? Or is calculating the normalizing constant itself your end objective? $\endgroup$ – jbowman Dec 5 '12 at 18:22
  • $\begingroup$ Thanks for your replies. @whuber I was coming to a similar conclusion, but your example clarifies that. In early iterations of the algorithm $f$ might be expected to have less variation - would MCMC be viable in that sort of situation? And if so, what sort of method could be used, as I'm having trouble finding eg Metropolis-Hastings over a discrete distribution? jbowman - Yes, the end objective is to choose an $x_i$ with weighted probability/sample from $f$, preferably without calculating the normalisation constant by evaluating $f$ at each $x_i$. $\endgroup$ – user17513 Dec 6 '12 at 9:22
2
$\begingroup$

Yes, you can generate random samples from $\{x_i\}_{i=1}^n$ where the associated probabilities $f(x_i)$ are known up to a multiplicative constant (the constant of integration (COI) using MCMC. However, you do have to be able to calculate $f(x_i)$ (without the COI) for all $i$ should $x_i$ be proposed. If your sample is large enough relative to $n$, you won't be saving yourself anything by not just calculating the COI by summing $f$ over all $i$.

Here's an example in R that does this for $n = 10$, using an independence sampler for the proposal (which simplifies the calculation of the acceptance probability $\alpha$). The proposal is for an index $j \in \{1, \dots, n\}$, although what is stored as the sampled value is the value of $x$ associated with the current index $i$ (if the proposal is rejected) or the accepted proposal $j$.

# Example x, fx
x <- c("a","b","c","d","e","f","g","h","i","j")
fx <- 0.01*(2^(0:9))   # Note this is unnormalized!

# i is the "current" index
# j is the "proposal" index
# y is the output stream of samples from x

# Initialization / allocation
y <- rep(0, 10000)
i <- sample(1:10,1)

# MCMC (Metropolis-Hastings)
u01 <- runif(length(y))
for (k in 1:length(y)) {
   j <- sample(1:10,1)  # Independence sampler
   alpha <- fx[j] / fx[i]
   if (u01[k] <= alpha) {
      y[k] <- x[j]
      i <- j
   } else {
      y[k] <- x[i]
   }
}

# Compare sample frequencies to true probabilities
SampleFreq <- table(y)/length(y)
df <- data.frame(Value=names(SampleFreq), 
                 Sample.Freq=as.numeric(SampleFreq), 
                 True.Prob=fx/sum(fx))
> df
   Value Sample.Freq True.Prob
1      a      0.0009 0.0009775
2      b      0.0022 0.0019550
3      c      0.0040 0.0039101
4      d      0.0093 0.0078201
5      e      0.0166 0.0156403
6      f      0.0345 0.0312805
7      g      0.0678 0.0625611
8      h      0.1282 0.1251222
9      i      0.2443 0.2502444
10     j      0.4922 0.5004888
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy