Closed Form for n-th digit Generalization in Benford’s Law

The sum as described in the Wikipedia article has a closed form. It took me some time but after some twiddling with the sum itself I decided to exponentiate the whole thing to get rid of the logarithm and had a simple product I was able to handle.

With d the digit and b the base:

{\displaystyle{ \prod_{n=b^{k-2}}^{b^{k-1}-1}\left(1+\frac{1}{d+b n}\right)  }}

This can be expressed as a fraction of two rising factorials (Pochhammer symbol).

{\displaystyle{ \frac{\Bigl(\frac{b^k+d b+b}{b^2}\Bigr)_{b^{k-2}\left( b - 1\right)}}{\Bigl(\frac{b^k+d b}{b^2}\Bigr)_{b^{k-2}\left( b - 1\right) }} }}

The rising factorial (x)_n can be expressed as a fraction of two Gamma functions

{\displaystyle{ (x)_n = \frac{\Gamma(x + n)}{\Gamma(x)}  }}

With which we get

{\displaystyle{\frac{\Gamma\bigl(b^{k-2}+\frac{d}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}+\frac{1}{b}\bigr)}{\Gamma\bigl(b^{k-2}+\frac{d}{b}+\frac{1}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}\bigr)}  }}

By setting d = b-1 we can already see that the result will get smaller the larger b gets with the obvious limit

{\displaystyle{ \lim_{b\to\infty} \frac{\Gamma\bigl(b^{k-2}+\frac{d}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}+\frac{1}{b}\bigr)}{\Gamma\bigl(b^{k-2}+\frac{d}{b}+\frac{1}{b}\bigr) \Gamma\bigl(b^{k-1}+\frac{d}{b}\bigr)} = 0 \quad\text{with }d = b-1 \land d\in\mathbb{N} }}

The arguments of the Gamma functions get quite large which makes the results hard to handle but we already need the logarithm of the product, that means we can use the logarithm of the Gamma function and save one step (we are on the real-line, so I used the equality \log\Gamma(x) = \log\left(\Gamma\left(x\right)\right) to avoid an even denser parentheses fence than it already is)

{\displaystyle{ \Biggl(\log\Gamma\biggl(b^{k-2}+\frac{d}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}+\frac{1}{b}\biggr)\Biggr) - \Biggl(\log\Gamma\biggl(b^{k-2}+\frac{d}{b}+\frac{1}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}\biggr)\Biggr) }}

We can simplify a bit with \frac{d}{b} + \frac{1}{b} = 1 because d = b-1. Furthermore b^{k-2} + \frac{d}{b} = \frac{b^{k-1} + d}{b} which can save us a multiplication here and there. So far we got

{\displaystyle{ \Biggl(\log\Gamma\biggl(\frac{b^{k-1}+d}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+1\biggr)\Biggr) - \Biggl(\log\Gamma\biggl(b^{k-2}+1\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}\biggr)\Biggr) }}


{\displaystyle{ \Biggl(\log\Gamma\biggl(\frac{b^{k-1}+d}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+1\biggr)\Biggr) - \Biggl(\log\Gamma\biggl(\frac{b^{k-1}+d + 1}{b}\biggr) + \log\Gamma\biggl(b^{k-1}+\frac{d}{b}\biggr)\Biggr) }}

In the special case I used it for b was a power of two such that none of the exponentiations had to be calculated, the numbers could be build in O(1).

Calling the whole thing f(d,b) we can compute the probability that the nth digit is b-1 with \frac{f(d,b)}{ \log(b)}. For general 0 \le d < b use the unsimplified version above for f(d,b).

Beware: the numbers get really large, you will need a slightly higher amount of precision for the loggamma function if you calculate with a large base. I had 226 as the base and had trouble to compute the values for places higher than the 20th in PARI/GP (it said *** lngamma: division by zero). It worked in calc with a precision of 300 decimal digits but it is sloooow. Who wrote that lame gamma function? Who? Me? Myself? Really? Ouch! šŸ˜‰


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s