Issue of precision when coding discrete distribution in c

Question

I am relatively new to C and so when it comes to issues of precision I am always a bit worried. I am working with a rather larger code and at a specific point I need to select a value based off of a discrete probability distribution. Essentially what I need to do is select a value k=1,2,...N with probabilities given by the formula

P(k) = [e^a*F(k)]/[Σ_i=1,...,Ne^aF(i)],

where F is just some function that uses k as an input and a is a real number. [Sorry the formula looks rough I wasn't sure how to format it.] My code looks (roughly) like this so far. I've had to make some changes to aid readability out of context.

int k;
int N=10;
double a=5.5; //In the code it is just a positive real number
double xtemp=0;
double X[10];
double fk;

for(k=0;k<N;k++){
  fk = F(k); // Here F returns a double
  xtemp += exp(a*fk);
  X[k]=xtemp;
}
xtemp = xtemp*U01(); //U01 returns a uniformly random number on [0,1]
for(k=0;k<K;k++){
   if(xtemp<=X[k]){
     return k;
   }
}

This should give the desired result from a mathematical perspective, however someone better at coding than me said the xtemp could become problematic and become imprecise. I was wondering how you could do error handling for this? I'd also be open to a better way to implement the idea. I appreciate any help and advice.

What are the individual `exp(a*F(k))` values like? Are they very big? Very small? Are they a mixture of very big and very small? The last one is most likely to cause precision issues ... — Useless, Jul 12 '17 at 16:59
This might help you: https://www.tutorialspoint.com/cprogramming/c_data_types.htm — Leonardo Alves Machado, Jul 12 '17 at 17:03
@Useless The values vary a good bit, and can end up being rather large. I'd say that typically 1 — Scott, Jul 12 '17 at 17:17
@Scot Have you tried to see whether `exp(720)` is even different from `INFINITY` in your C implementation? — EOF, Jul 12 '17 at 17:18
`exp(720)` is way too big for a `double` but apparently `long double` should be able to handle it. — Patrick87, Jul 12 '17 at 17:36
@EOF I have tried it and it registers `exp(720)` as infinity. — Scott, Jul 12 '17 at 17:38
@Scot I'd say that confirms that the code you have provided is *not* acceptable. — EOF, Jul 12 '17 at 17:39
Unsigned doubles are not a thing, btw, and one bit wouldn't make that much of a difference anyway — Useless, Jul 12 '17 at 18:11
Oh, and since you're trying to store the cdf values in your array, have you tried just calculating this in python (or something else with bignums) ahead of time, and seeing if you can scale the values back down while keeping the proportions? — Useless, Jul 12 '17 at 18:13
@Patrick87 There is no guarantee `long double` can represent it. — too honest for this site, Jul 12 '17 at 19:48

score 0 · Answer 1 · answered Jul 12 '17 at 17:48

Given that your current code could cause evaluation of sub-expressions that fall outside of the finite range of representable values for double, I'd recommend rewriting the expression. Consider something like this:

P(k) = [e^a*F(k)]/[Σ_i=1,...,Ne^{a * F(i)}]

P^-1(k) = [Σ_i=1,...,Ne^{a * F(i)}] / [e^a*F(k)]

P^-1(k) = Σ_i=1,...,N[e^{a * F(i)} / e^a*F(k)]

P^-1(k) = Σ_i=1,...,N[e^{a * F(i) - a * F(k)}]

P^-1(k) = Σ_i=1,...,N[e^{a * (F(i) - F(k))}]

P(k) = [Σ_i=1,...,N[e^{a * (F(i) - F(k))}]]^-1

If you could provide a minimal, complete and testable example, you could probably get considerably better advice.

Issue of precision when coding discrete distribution in c

1 Answers1