I have a pdf that I need to sample from, but I'm having numerical issues and was hoping someone could provide some insight. My pdf is a product of gaussians

$$ P(\theta |d ) = \prod_{i=1}^{N} C_{i} e^{-\frac{1}{2 \sigma_{i}^{2}}(d_{i} - y_{i}(\theta) )^{2} } \,. $$

Lets say for instance that $\theta$ is only one dimensional, making $P(\theta | d)$ a one dimensional pdf. I wish to plot P(\theta | d) as a function of $\theta$, where $\theta$ is on the x-axis of the plot and $P(\theta | d)$ is on the y-axis.

In this case I can simply plug in a bunch of various values for $\theta$ and plot the result. The problem is, the value inside the exponent which I'll call $K$,

$$ K = - \sum_{i = 1}^{N} \frac{1}{2 \sigma^{2}_{i}}(d_{i} - y_{i}(\theta) )^{2} $$

can vary very much, and never really approaches zero. For example, the range of values for K I got varies from about -20000 to about -500. This is because my data, $d_{i}$, is pretty fairly well corrupted by noise.

Is there anything I can do to plot this pdf as a function of $\theta$? Obviously everything I do on my computer (taking $e^{K}$ and getting mad) is pretty much giving me all zeros.


If your goal is to look at the data and get a feel for what's going on in the plot, you may be stuck. Exponents varying by that much I don't think are going to plot nicely. Your best bet might be to plot $K$ instead, as it's at least close to being the log probability.

If your goal is to draw samples though, have you considered using metropolis hastings? I ask, as the only part of that which interacts with the distribution in question (i.e. the part where the algorithm decides whether to accept a new sample $x'$ given the prior sample $x$) depends on the ratio of probabilites $a$, where $a = P(x')/P(x)$.

But to get the ratio, you don't need to calculate the actual probabilities $P(x')$ and $P(x)$. Instead you can use

$$ a = \frac{P(x')}{P(x)} = e^{\log(P(x')) - \log(P(x))}. $$

That means you can just calculate the log probabilities directly (which should be scaled much more nicely), subtract them (getting rid of any huge nasty constant factors between the two), and exp the answer. You won't be free of numeric issues, but they ought to be lessened.

  • $\begingroup$ Clever! Great answer. All I could think of is trying to jam this to fit the logsumexp trick, but this is just right. $\endgroup$
    – bill_e
    Feb 28 '14 at 22:59

plot K itself instead of $e^K$, or plot $\log_{10}(e^K)$ if you want the plot be easier to read.

  • $\begingroup$ OP can't literally take $log_{10} e^K$ since $e^K$ numerically underflows. $\endgroup$
    – Dave
    Feb 28 '14 at 16:55
  • 1
    $\begingroup$ $\log_{10}e^K=K\log_{10}e$, it's basically the same as K, but much easier to interpret since it's in powers of 10 $\endgroup$
    – Aksakal
    Feb 28 '14 at 16:57
  • $\begingroup$ No, this isn't suitable for me. Of course there are many ways to not plot the pdf. $\endgroup$
    – bill_e
    Feb 28 '14 at 22:57
  • $\begingroup$ it's a plot PDF, but in Log scale. $\endgroup$
    – Aksakal
    Feb 28 '14 at 23:10

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.