Home / Previous Lesson / Next Lesson
Related Dyspoissometer macros: LOG(), LOG_N_PLUS_1(), LOG_SUM(), LOG_SUM_N_PLUS_1().
UPDATE: Interval math considerations, and a suggested accelerated approximation method, are discussed in this whitepaper which was also published here at vixra.org. sha256sum=39895297e3b50093bcc6b0595f2e1d1cf2639a8ffa7bb942487c0256dbff5b8f
NOTE: If you arrived here because you searched for 0.918938533204672741780..., then you are dealing with ((1/2) log(2π)), which suggests that your math problem may have something to do with logarithmic summation.
This lesson is about how to compute the natural logarithm ("log") of a factorial, namely;
log(N!)
where N is an integer greater than 1. (In the trivial cases, where N is zero or one, then (N!) is one by definition, so L is zero.) We refer to log(N!) as the "logsum" of N (which is the sum of the logs, not the log of the sum).
Clearly, we can evaluate it by brute force:
log(N!) = Σ(M = (1, N), log(M))
This works quickly and accurately for small N, but at a certain point, neither adjective applies. Fortunately, we can use what I call "Bayes' Approximation" to provide a huge extension in the domain of N which satisfies both adjectives.
A bit of background: James Stirling popularized an approximation for the logsum which was first published by Abraham de Moivre sometime prior to his death in 1754. This became known as "Stirling's Approximation", even though it was neither his discovery nor, as it turns out, convergent!
To make a long story short, Stirling's Approximation comes from approximating log(N!) as an integral using calculus, instead of what it literally is, which is the sum of (a finite number of (an infinite number of rational fractions)). (The log itself is computable using various infinite series, the most rapidly converging of which being the "area hyperbolic tangent" series.) Stirling (well, de Moivre) used the fact that:
∫ log(x) dx = x log(x) - x
if we take the constant of indefinite integration as zero. Integrating this function on the domain [1, N] produces an approximation of log(N!). Then, by adding evermore corrective terms, the approximation can approach the exact value.
At least, that was the theory. But in 1763, 2 years after Thomas Bayes' death, the Royal Society published a letter from him to John Canton, in which he demonstrated that the approximation which had been misappropriated to Stirling was in any event nonconvergent! In other words, after a certain number of terms, it became increasingly less accurate. Oops.
I don't know who rectified the situation. Perhaps it was Stirling himself, but in my view, Bayes' deserves the credit for finding the problem, and certainly we want to draw a distinction between "Stirling's Broken Approximation" and "Someone's Fixed Approximation". It is the latter which I call "Bayes' Approximation". Yet for some vexing reason, Stirling's Broken Approximation appears to be in popular use today, perhaps because the first few terms are similar to those of Bayes' (be careful!). I wonder how many products have been misengineered as a result...
Now it turns out that Bayes' Approximation, which is to my knowledge absolutely convergent for (N > 0), has buried within it a deep combinatorics problem not unlike the determination of way count itself! I suppose this logical, since logfreedom, which contains logsum terms, is informationally equivalent to way count. Here is how the approximation works:
log(N!) = (N + (1/2)) log(N + 1) - (N + 1) + (1/2) log(2π) + C(N)
All of the above is derived from calculus with an optimized constant of indefinite integration. Note that the approximation is officially for computing log(Γ(N)), which I have adapted to logsum using the following identity:
log(N!) ≡ log(Γ(N+1))
For its part, the C version of Dyspoissometer uses lgamma(N + 1) to compute the log(N!) in double precision and lgammaq() to compute it in quad precision; the Mathematica version uses LogGamma[] (as opposed to Log[Gamma[]], which is slow).
But what about C(N)? It's given by the following infinite series of finite products:
C(N) = Σ(M = (1, ∞), c(M) / Π(P = (1, M), (N + P + 1)))
where c(M) is the finite series given by:
c(M) = (1 / (2M)) Σ(P = (1, M), P |s(M, P)| / ((P + 1)(P + 2)))
where |s(M, P)| (note the bars) is the absolute value of the Stirling numbers of the first kind (which, strictly speaking, are signed). Those numbers are easily generated, if you trudge through Wikipedia's awkward presentation of the algo. Without attempting to duplicate the article here, we can see the echos of way count in their mathematical properties: "In particular, the Stirling numbers of the first kind count permutations according to their number of cycles (counting fixed points as cycles of length one)." So perhaps there is a connection from way count to graph theory which I have not yet appreciated.
Putting it all together:
log(N!) = (N + (1/2)) log(N + 1) - (N + 1) + (1/2) log(2π) + (1 / (12 (N + 2))) + (1 / (12 (N + 2) (N + 3))) + (59 / (360 (N + 2) (N + 3) (N + 4))) + (29 / (60 (N + 2) (N + 3) (N + 4) (N + 5))) + ...
Quiz:
Find log(7!) using brute force and Bayes' Approximation.
Answer:
1. log(7) + log(6) + log(5) + log(4) + log(3) + log(2) = 8.525...
2. Using Bayes' Approximation, we get approximately the same answer at WolframAlpha here.
Fortunately -- at least as far as I can tell -- the lgamma() function in GCC uses Bayes' Approximation. (I discovered this after painstakingly implementing the above formula for a week, only to discover that lgamma() was both slightly slower and more accurate, yet radically simpler by virtue of being a builtin function. All this was quite a coincidence, triggered by opening to the wrong instruction manual page, whereupon I discovered lgamma.)
By the way, there other approaches to approximating the logsum, none of which I prefer to Bayes' Approximation. But see for yourself here at StackExchange.
Now that you know how to compute the logsum efficiently, we can explore the applications of logsum in the next lesson.
Keywords: 0.918938533204672741780, 0.91893853320467274178, 0.9189385332046727417, 0.918938533204672741, 0.91893853320467274, 0.9189385332046727, 0.918938533204672, 0.91893853320467, 0.9189385332046, 0.918938533204, 0.91893853320, 0.9189385332, 0.918938533, 0.91893853, 0.9189385, 0.918938, 0.91893
Related Dyspoissometer macros: LOG(), LOG_N_PLUS_1(), LOG_SUM(), LOG_SUM_N_PLUS_1().
UPDATE: Interval math considerations, and a suggested accelerated approximation method, are discussed in this whitepaper which was also published here at vixra.org. sha256sum=39895297e3b50093bcc6b0595f2e1d1cf2639a8ffa7bb942487c0256dbff5b8f
NOTE: If you arrived here because you searched for 0.918938533204672741780..., then you are dealing with ((1/2) log(2π)), which suggests that your math problem may have something to do with logarithmic summation.
This lesson is about how to compute the natural logarithm ("log") of a factorial, namely;
log(N!)
where N is an integer greater than 1. (In the trivial cases, where N is zero or one, then (N!) is one by definition, so L is zero.) We refer to log(N!) as the "logsum" of N (which is the sum of the logs, not the log of the sum).
Clearly, we can evaluate it by brute force:
log(N!) = Σ(M = (1, N), log(M))
This works quickly and accurately for small N, but at a certain point, neither adjective applies. Fortunately, we can use what I call "Bayes' Approximation" to provide a huge extension in the domain of N which satisfies both adjectives.
A bit of background: James Stirling popularized an approximation for the logsum which was first published by Abraham de Moivre sometime prior to his death in 1754. This became known as "Stirling's Approximation", even though it was neither his discovery nor, as it turns out, convergent!
To make a long story short, Stirling's Approximation comes from approximating log(N!) as an integral using calculus, instead of what it literally is, which is the sum of (a finite number of (an infinite number of rational fractions)). (The log itself is computable using various infinite series, the most rapidly converging of which being the "area hyperbolic tangent" series.) Stirling (well, de Moivre) used the fact that:
∫ log(x) dx = x log(x) - x
if we take the constant of indefinite integration as zero. Integrating this function on the domain [1, N] produces an approximation of log(N!). Then, by adding evermore corrective terms, the approximation can approach the exact value.
At least, that was the theory. But in 1763, 2 years after Thomas Bayes' death, the Royal Society published a letter from him to John Canton, in which he demonstrated that the approximation which had been misappropriated to Stirling was in any event nonconvergent! In other words, after a certain number of terms, it became increasingly less accurate. Oops.
I don't know who rectified the situation. Perhaps it was Stirling himself, but in my view, Bayes' deserves the credit for finding the problem, and certainly we want to draw a distinction between "Stirling's Broken Approximation" and "Someone's Fixed Approximation". It is the latter which I call "Bayes' Approximation". Yet for some vexing reason, Stirling's Broken Approximation appears to be in popular use today, perhaps because the first few terms are similar to those of Bayes' (be careful!). I wonder how many products have been misengineered as a result...
Now it turns out that Bayes' Approximation, which is to my knowledge absolutely convergent for (N > 0), has buried within it a deep combinatorics problem not unlike the determination of way count itself! I suppose this logical, since logfreedom, which contains logsum terms, is informationally equivalent to way count. Here is how the approximation works:
log(N!) = (N + (1/2)) log(N + 1) - (N + 1) + (1/2) log(2π) + C(N)
All of the above is derived from calculus with an optimized constant of indefinite integration. Note that the approximation is officially for computing log(Γ(N)), which I have adapted to logsum using the following identity:
log(N!) ≡ log(Γ(N+1))
For its part, the C version of Dyspoissometer uses lgamma(N + 1) to compute the log(N!) in double precision and lgammaq() to compute it in quad precision; the Mathematica version uses LogGamma[] (as opposed to Log[Gamma[]], which is slow).
But what about C(N)? It's given by the following infinite series of finite products:
C(N) = Σ(M = (1, ∞), c(M) / Π(P = (1, M), (N + P + 1)))
where c(M) is the finite series given by:
c(M) = (1 / (2M)) Σ(P = (1, M), P |s(M, P)| / ((P + 1)(P + 2)))
where |s(M, P)| (note the bars) is the absolute value of the Stirling numbers of the first kind (which, strictly speaking, are signed). Those numbers are easily generated, if you trudge through Wikipedia's awkward presentation of the algo. Without attempting to duplicate the article here, we can see the echos of way count in their mathematical properties: "In particular, the Stirling numbers of the first kind count permutations according to their number of cycles (counting fixed points as cycles of length one)." So perhaps there is a connection from way count to graph theory which I have not yet appreciated.
Putting it all together:
log(N!) = (N + (1/2)) log(N + 1) - (N + 1) + (1/2) log(2π) + (1 / (12 (N + 2))) + (1 / (12 (N + 2) (N + 3))) + (59 / (360 (N + 2) (N + 3) (N + 4))) + (29 / (60 (N + 2) (N + 3) (N + 4) (N + 5))) + ...
Quiz:
Find log(7!) using brute force and Bayes' Approximation.
Answer:
1. log(7) + log(6) + log(5) + log(4) + log(3) + log(2) = 8.525...
2. Using Bayes' Approximation, we get approximately the same answer at WolframAlpha here.
Fortunately -- at least as far as I can tell -- the lgamma() function in GCC uses Bayes' Approximation. (I discovered this after painstakingly implementing the above formula for a week, only to discover that lgamma() was both slightly slower and more accurate, yet radically simpler by virtue of being a builtin function. All this was quite a coincidence, triggered by opening to the wrong instruction manual page, whereupon I discovered lgamma.)
By the way, there other approaches to approximating the logsum, none of which I prefer to Bayes' Approximation. But see for yourself here at StackExchange.
Now that you know how to compute the logsum efficiently, we can explore the applications of logsum in the next lesson.
Keywords: 0.918938533204672741780, 0.91893853320467274178, 0.9189385332046727417, 0.918938533204672741, 0.91893853320467274, 0.9189385332046727, 0.918938533204672, 0.91893853320467, 0.9189385332046, 0.918938533204, 0.91893853320, 0.9189385332, 0.918938533, 0.91893853, 0.9189385, 0.918938, 0.91893
No comments:
Post a Comment