CODEWITHSUNDEEP
python-3.xfloating-pointprecisionbayesianprobability-distribution
Aditya Agarwal
Aditya Agarwal

Reputation: 231

Calculating a custom probability distribution in python (numerically)

I have a custom (discrete) probability distribution defined somewhat in the form: f(x)/(sum(f(x')) for x' in a given discrete set X). Also, 0<=x<=1. So I have been trying to implement it in python 3.8.2, and the problem is that the numerator and denominator both come out to be really small and python's floating point representation just takes them as 0.0.
After calculating these probabilities, I need to sample a random element from an array, whose each index may be selected with the corresponding probability in the distribution. So if my distribution is [p1,p2,p3,p4], and my array is [a1,a2,a3,a4], then probability of selecting a2 is p2 and so on.
So how can I implement this in an elegant and efficient way?
Is there any way I could use the np.random.beta() in this case? Since the difference between the beta distribution and my actual distribution is only that the normalization constant differs and the domain is restricted to a few points.

Note: The Probability Mass function defined above is actually in the form given by the Bayes theorem and f(x)=x^s*(1-x)^f, where s and f are fixed numbers for a given iteration. So the exact problem is that, when s or f become really large, this thing goes to 0.

Upvotes: 0

Views: 206

Answers (1)

dmuir
dmuir

Reputation: 4431

You could well compute things by working with logs. The point is that while both the numerator and denominator might underflow to 0, their logs won't unless your numbers are really astonishingly small.

You say

f(x) = x^s*(1-x)^t

so

logf (x) = s*log(x) + t*log(1-x)

and you want to compute, say

p = f(x) / Sum{ y in X | f(y)}

so

p = exp( logf(x) - log sum { y in X | f(y)}
  = exp( logf(x) - log sum { y in X | exp( logf( y))}

The only difficulty is in computing the second term, but this is a common problem, for example here

On the other hand computing logsumexp is easy enough to to by hand.

We want

S = log( sum{ i | exp(l[i])})

if L is the maximum of the l[i] then

S = log( exp(L)*sum{ i | exp(l[i]-L)})
  = L + log( sum{ i | exp( l[i]-L)})

The last sum can be computed as written, because each term is now between 0 and 1 so there is no danger of overflow, and one of the terms (the one for which l[i]==L) is 1, and so if other terms underflow, that is harmless.

This may however lose a little accuracy. A refinement would be to recognize the set A of indices where

l[i]>=L-eps (eps a user set parameter, eg 1)

And then compute

N = Sum{ i in A | exp(l[i]-L)}
B = log1p( Sum{ i not in A | exp(l[i]-L)}/N)
S = L + log( N) + B

Upvotes: 2

Related Questions