# HP Forums

Full Version: [WP34s] Inverse discrete probability distributions
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

Hello all,

The post will be long, but I see no other way to make what I think is a very important point clearly enough.

I was helping my daughter with her homework. She had to spin a roulette many times and count the frequencies of different outcomes. I wrote a 34s program to do a Monte Carlo simulation. That way she could do about 1000 rolls in 2 seconds (a loop calling RAN# and increasing the appropriate register each time).

Just for fun, I decided to write a second program to do the same but using the inverse binomial distribution instead, only to quickly realize that the implementation of inverse discrete distributions on the 34s is "off by one".

Let me explain. For simplicity let us consider the geometric distribution, say with parameter p=0.5 stored in J.

Such a geometric variable takes the value m=0 with probability 1/2, m=1 with probability 1/4, m=2 with prob=1/8, etc.

The distribution function (the function Geom in the 34s), is the function F(x) defined for all real values of x as follows:

F(x) = 0 if x<0

F(x) = 1/2 if 0<=x<1

F(x) = 3/4 if 1<=x<2

F(x) = 7/8 if 2<=x<3

...

F(x) = 1 - 1/2^n if n-1<=x<n

...

and so on ad infinitum. It's OK to put F(-infinity)=0 and F(+infinity)=1 if one is so inclined.

The "inverse distribution" F^(-1)(x), that is, the inverse function to F(x) above is, strictly speaking only defined for all real values of x, but only for the values of x=0,1/2,3/4,7/8,..., etc. And worse, it is not *uniquely* defined for those values. For instance, we could take F^(-1)(1/2) to be equal to any number x such that 0<=x<1. This is mathematically correct, strictly speaking.

The question becomes: What is the more *useful* way to define the inverse functions?

The 34s calls that inverse distribution function above Geom-1. Currently, the 34s calculates this function as follows:

Geom-1(x) = 0 if x<0,

Geom-1(x) = -1 if 0<=x<1/2,

Geom-1(x) = 0 if 1/2<=x<3/4,

Geom-1(x) = 1 if 3/4<=x<7/8,

Geom-1(x) = 2 if 7/8<=x<15/16,

etc.

!!!

The output for x<0 notwithstanding (which I propose to make simply undefined), the values returned for x>=0 are almost always one unit smaller than they ought to be. The values returned should be:

Geom-1(x) is undefined if x<0 (one can probably argue this should be either -1, or -infinity, or NaN. The point is that nobody should ever need to evaluate an inverse distribution function when the argument is not in the range 0<=x<=1.),

Geom-1(x) = 0 if 0<=x<=1/2,

Geom-1(x) = 1 if 1/2<x<=3/4,

Geom-1(x) = 2 if 3/4<x<=7/8,

Geom-1(x) = 3 if 7/8<x<=15/16,

etc.

Unfortunately, it seems as if all of the discrete inverse distributions are "off by 1" in the 43s currently (Caveat: I only checked the geometric and binomial).

I want to point out an important immediate corollary of fixing the inverse distributions as I propose above: One can use the inverse distribution (that is, the "inverse transform method" as explained in these notes) to go from a uniform RAN# variable to a variable with the desired distribution:

RAN#

Geom-1

generates a random variable distributed geometrically. This method fails with the current implementation of the inverse discrete random variables. If we replace Geom-1 by, say, Binom-1, Expon-1 or Phi-1 we get a binomially, exponentially or normally distributed random variable, regardless of whether the distribution is discrete or continuous. Elegant and consistent.

Eduardo

PS: I have the source code of the probability distributions, but I could not quickly make sense of it so as to point out more specifically how to implement the changes. Sorry.

[Edited to exchange some <= to < in order that the examples I give indeed be inverse to the original distribution.]

Edited: 15 May 2012, 11:40 p.m. after one or more responses were posted

Eduardo,

Thanks for the link and for reporting your observations. It will take us some time to check and test, but we will return to you as soon as we know more :-)

Actually, the link I provided is probably not a very good reference on the "inverse distribution method", but merely the first document a quick search turned up.

However, this is a standard result found in just about any textbook on prob/stats.

If you want me to send a more detailed/technical explanation (or a better reference) I can provide it, though perhaps the forum may not be the best place for the discussion of the details.

(Can always reach me by e-mail: firstname dot lastname at gmail dot com)

Eduardo

This is by design. The inverse of the discrete functions returns, for a given probability p and distribution X, the largest integer n such that the Pr(X <= n) <= p. I.e. it is rounding the value associated with the probability down. Thus for 0 which has Pr(X=0) of 0.5, the inverse won't return 0 until p >= 0.5.

What you seem to want is for the inverse to round up. I.e. the smallest integer n such that Pr(x <= n) > p.

Both approaches are valid, the question is which one is more desirable here. Anyone want to make a case either way?

The x<0 case is a bug and it will be fixed once this question is answered.

- Pauli

Edited: 15 May 2012, 9:03 p.m.

Quote:
PS: I have the source code of the probability distributions, but I could not quickly make sense of it so as to point out more specifically how to implement the changes. Sorry.

The final fix for the discrete distributions is at label qf_discrete_final in trunk/xrom/distributions/geometric.wp34s.

Changing this is the easiest way to alter the rounding behaviour of all discrete distributions at once.

The question remains: what behaviour is most desirable?

- Pauli

Well, my original post advocates the "rounding up" method, and my opinion is strongly for it.

Although I'm no Wolfram fanboy, they have a decent discussion on the "method of inverse transforms" http://demonstrations.wolfram.com/TheMethodOfInverseTransforms/

The "Details" section of that page explains how to use the inverse distribution function to generate a variable with any given distribution starting from a uniformly distributed variable as simulated by RAN#. I want to emphasize that this method *only* works if the inverse distribution is computed rounding up rather than down.

I consider the use of the inverse distribution function in transforming uniformly distributed RAN# numbers into random numbers with a *different* distribution to be the quintessential use of these inverses (when "rounding up" is used, as I have advocated, in the way I find mathematically more natural). So I ask the question: Under what circumstances would it be more natural or useful to define the inverse distribution by rounding down?

If I have a chance tomorrow, I'll try to find a proper reference (e.g., a textbook) on this matter.

Eduardo

For reference, here's a couple of web pages describing how the inverse distribution function should be defined (either of these is far better than the reference I gave in my original post, which although correct is just lecture notes for a sophomore-level course in actuarial models).

First one, by faculty at the Computer Science Dept. of Ben Gurion University in Israel:

Another one, written in complete generality (using slightly higher mathematical lingo and including a proof of the Inverse Distribution Theorem), by Ursa Pantle, a statistician at Universität Ulm in Germany:

In a nutshell, for a variable with distribution F(x) = probability that the random variable takes values <=x, the generalized inverse distribution function is defined by:

F-1(y) = infimum of values x such that F(x)>=y.

In particular, if F(x) < y then F-1(y) must be at least x, and in fact if the underlying distribution is discrete, then F-1(y)>x ---in other words, rounding must be done up, never down.

Eduardo

These links are talking about generating distributed random numbers -- this doesn't necessary relate to the inverse cumulative distribution function for a discrete random variable. The inverse function is ill defined -- both what the 34S and what you want are inverses in the sense that they return the correct values over the integers which is the domain of the discrete distributions.

That all said, nobody else has said anything thus far and I really don't mind which way these functions are defined (it is only an INC or DEC between them), so I'll go ahead and change them to the upward rounding flavours in due course.

- Pauli

Quote:
That all said, nobody else has said anything thus far and I really don't mind which way these functions are defined (it is only an INC or DEC between them), so I'll go ahead and change them to the upward rounding flavours in due course.

Precisely, since the strict inverses are undefined for most real numbers, it becomes a matter of extending the definition in the most useful manner. Probabilists and statisticians have a preferred way of defining this generalized inverse (viz., the documents I have referenced in my prior posts in this thread). I noticed the 34s inverse discrete distribution functions differ from such definition in most of their domain, and I have tried pointing out so. That said, I certainly do not wish to impose the work on you or the 34s team. I'm sure I can figure out how to make the code changes myself and contribute them to the project given a little time. So long as the changes are welcome, that is.

Eduardo

The changes are already done, just not built yet. They'll need testing however...

The only constraint I remember on discrete distributions was they had to be right continuous and both of these definitions are that. Well, beyond always being non-negative and integrating/summing to unity as all distributions must. It's been a while since I classed myself a statistician and longer still since I studied postgraduate probability theory so I might have forgotten something....

Anyway, the point is moot now.

- Pauli

Quote:
The changes are already done, just not built yet. They'll need testing however...

Consider me first in line to test! (In due course.)

Eduardo

Quote:
Consider me first in line to test! (In due course.)

Well, you can test! ;-)