Tuesday, November 3, 2015

Stirling's Approximation: Derivation and Corollaries

 

Lemma 1: \(\lim_{n \rightarrow \infty} \sqrt[n]{n!}/n=1/e\)


Way 1 (somewhat rigorous)


From elementary calculus, we have that: \[ \int_{0}^{1} \ln(x) dx =-1 \] Taking this as a Riemann sum, as done in introductory calculus, we have: \[ -1=\int_{0}^{1}\ln(x)dx=\lim_{N \rightarrow \infty} \sum_{k=1}^{N}\ln\left (\frac{k}{N} \right ) \cdot \frac{1}{N} \] \[ -1=\lim_{N \rightarrow \infty} -\ln(N)+\frac{1}{N} \sum_{k=1}^{N}\ln\left (k \right ) \] \[ -1=\lim_{N \rightarrow \infty} -\ln(N)+\frac{1}{N} \ln\left (N! \right ) \] Therefore, \[ \lim_{N \rightarrow \infty} \frac{\sqrt[N]{N!}}{N}=\frac{1}{e} \]



Way 2 (less rigorous)


\[ \lim_{n \rightarrow \infty} \frac{\sqrt[n]{n!}}{n}=x \] So, for n big, in a certain sense: \[ n! \approx (nx)^n \] \[ \frac{(n+1)!}{n!(n+1)}=1 \approx \frac{((n+1)x)^{n+1}}{(nx)^n (n+1)}=\left ( 1+ \frac{1}{n} \right )^n x \] Thus, in order to get equality in the limit, we must have: \[ x = \lim_{n \rightarrow \infty} \left ( 1+ \frac{1}{n} \right )^{-n}=\frac{1}{e} \]



Lemma 2: Wallis Product in Factorial Form


Recall from this article the following expression for pi: \[ \frac{\pi}{2}=\prod_{k=1}^{\infty}\frac{2k \cdot 2k}{(2k-1)(2k+1)}=\lim_{N \rightarrow \infty}\prod_{k=1}^{N}\frac{2k \cdot 2k}{(2k-1)(2k+1)}=\lim_{N \rightarrow \infty} \frac{\left ( 2^N \cdot N! \right )^4}{\left ( (2N)! \right )^2(2N+1)} \]



Lemma 3: An Inequality for the Natural Logarithm


Let \(x,y > 0\). Clearly \[ 0 \leq \frac{1}{y^2 (1+y)^2 (2y+1)^2} \] Therefore \[ 0 \leq \int_{x}^{\infty}\frac{dy}{y^2 (1+y)^2 (2y+1)^2}=\frac{1}{x}+\frac{1}{x+1}+\frac{4}{x+1/2}-6\ln \left ( 1+\frac{1}{x} \right ) \] \[ 6\ln \left ( 1+\frac{1}{x} \right ) -\frac{6}{x+1/2} \leq \frac{1}{x}+\frac{1}{x+1}-\frac{2}{x+1/2} \] \[ (x+\tfrac{1}{2})\ln \left ( 1+\frac{1}{x} \right ) -1 \leq \frac{(x+\tfrac{1}{2})}{6}\left (\frac{1}{x}+\frac{1}{x+1} \right )-\frac{1}{3}=\frac{1}{12x(x+1)} \] Also, clearly \[ 0 \leq \frac{16y^2+41y+24}{y(1+y)^2 (2+y)^2 (2y+1)^2} \] Therefore \[ 0 \leq \int_{x}^{\infty}\frac{16y^2+41y+24}{y(1+y)^2 (2+y)^2 (2y+1)^2} dy=6\left (\ln \left ( 1+\frac{1}{x} \right )-\frac{1}{x+\tfrac{1}{2}} \right)-\frac{1}{2(x+\tfrac{1}{2})(x+1)(x+2)} \] And so \[ \frac{1}{12(x+1)(x+2)} \leq (x+\tfrac{1}{2})\ln \left ( 1+\frac{1}{x} \right )-1 \]



Theorem: Stirling's Approximation


Let us define a function and sequence of coefficients as follows: \[ g(n)=\ln\left ( \frac{n!}{\left ( \tfrac{n}{e} \right )^n \sqrt{2\pi n}} \right )=\sum_{k=-\infty}^{\infty} A_k n^k \] We then have, from lemma 1, \[ \frac{1}{e}=\lim_{n \rightarrow \infty} \frac{\sqrt[n]{n!}}{n}=\lim_{n \rightarrow \infty} \frac{\sqrt[n]{\left (\tfrac{n}{e} \right )^n \sqrt{2\pi n} \cdot e^{g(n)}}}{n}=\frac{1}{e} \lim_{n \rightarrow \infty} \sqrt[2n]{2\pi n} \cdot e^{g(n)/n} \] Thus \[ 1=\lim_{n \rightarrow \infty} e^{g(n)/n}=\exp\left (\lim_{n \rightarrow \infty} \sum_{k=-\infty}^{\infty} A_k n^{k-1} \right )=\exp\left (\lim_{n \rightarrow \infty} \sum_{k=1}^{\infty} A_k n^{k-1} \right ) \] And therefore \(A_k=0\) for \(k \geq 1\). From lemma 2, \[ \frac{\pi}{2}=\lim_{n \rightarrow \infty} \frac{\left ( 2^n \cdot n! \right )^4}{\left ( (2n)! \right )^2(2n+1)}=\lim_{n \rightarrow \infty} \frac{\left ( 2^n \cdot \left (\tfrac{n}{e} \right )^n \sqrt{2\pi n} \cdot e^{g(n)} \right )^4}{\left ( \left (\tfrac{2n}{e} \right )^{2n} \sqrt{4\pi n} \cdot e^{g(2n)} \right )^2(2n+1)} \] \[ \frac{\pi}{2}=\lim_{n \rightarrow \infty} \frac{\left ( 2^n \cdot n! \right )^4}{\left ( (2n)! \right )^2(2n+1)}=\lim_{n \rightarrow \infty} \frac{\left ( 2^n \cdot \left (\tfrac{n}{e} \right )^n \sqrt{2\pi n} \cdot e^{g(n)} \right )^4}{\left ( \left (\tfrac{2n}{e} \right )^{2n} \sqrt{4\pi n} \cdot e^{g(2n)} \right )^2(2n+1)} \] \[ \frac{\pi}{2}=\lim_{n \rightarrow \infty} \frac{n \pi}{2n+1} \cdot e^{4g(n)-2g(2n)} \] \[ 0=\lim_{n \rightarrow \infty} 2g(n)-g(2n)=2A_0-A_0=A_0 \] Therefore, \(A_k=0\) for \(k \geq 0\), and thus \(\lim_{n \rightarrow \infty} g(n)=0\). Thus it follows that \[ \lim_{n \rightarrow \infty} \frac{n!}{\left ( \tfrac{n}{e} \right )^n \sqrt{2\pi n}}=1 \] This fact is known as Stirling's Approximation. Moreover, we have \[ g(n)-g(n+1)=\ln\left ( \frac{n!\left ( \tfrac{n+1}{e} \right )^{n+1} \sqrt{2\pi (n+1)}}{(n+1)!\left ( \tfrac{n}{e} \right )^n \sqrt{2\pi n}} \right )=\ln\left ( \frac{(n+1)^{n+\tfrac{1}{2}}}{e \cdot n^{n+\tfrac{1}{2}}} \right ) \] \[ g(n)-g(n+1)=(n+\tfrac{1}{2})\ln\left ( 1+\frac{1}{n} \right )-1 \] By lemma 3, we then have \[ \frac{1}{12(n+1)(n+2)} \leq g(n)-g(n+1) \leq \frac{1}{12n(n+1)} \] \[ \sum_{k=n}^{\infty} \frac{1}{12(k+1)(k+2)}=\frac{1}{12(n+1)} \leq \sum_{k=n}^{\infty} g(k)-g(k+1)=g(n)-g(\infty)=g(n) \leq \sum_{k=n}^{\infty} \frac{1}{12k(k+1)}=\frac{1}{12n} \] That is \(\tfrac{1}{12(n+1)} \leq g(n) \leq \tfrac{1}{12n}\). And therefore: \[ \left (\tfrac{n}{e} \right )^n \sqrt{2\pi n}\cdot e^{\tfrac{1}{12(n+1)}} \leq n! \leq \left (\tfrac{n}{e} \right )^n \sqrt{2\pi n} \cdot e^{\tfrac{1}{12n}} \] In fact, it is possible to obtain exact formulas for \(g(n)\). For example, by more advanced calculations, we can show that \[ g(n)=\int_{0}^{\infty}\frac{2 \tan^{-1}\left ( \tfrac{y}{n} \right )}{e^{2\pi y}-1}=\sum_{k=1}^{\infty} \frac{B_{2k}}{2k(2k-1)n^{2k-1}}=\frac{1}{12n}-\frac{1}{360n^3}+\frac{1}{1260n^5}- \cdots \] Where \(B_m\) is the mth Bernoulli number. These two expressions are, respectively Binet's second expression and Stirling's series.



Corollary: Product of a Rational Function


Firstly, since \[ \prod_{k=1}^N \left(ak+b\right)=a^N\prod_{k=1}^N \left(k+\frac{b}{a}\right) \] We will just evaluate \[ \prod_{k=1}^N \left(k+b\right)=\frac{(N+b)!}{b!} \approx \left(\frac{N+b}{e}\right)^{N+b}\frac{\sqrt{2\pi(N+b)}}{b!}=N^{N+b+\tfrac{1}{2}}e^{-N}\frac{\sqrt{2\pi}}{b!}e^{-b}\left(1+\frac{b}{N}\right)^{N+b+\tfrac{1}{2}} \] \[ \prod_{k=1}^N \left(k+b\right)=\frac{(N+b)!}{b!} \approx N^{N+b+\tfrac{1}{2}}e^{-N}\frac{\sqrt{2\pi}}{b!} \] More generally, given the above, it is not difficult to demonstrate the following generalization. Let \(m,n > 0\). Let \(a_1,a_2,...,a_m\) and \(b_1,b_2,...,b_n\) and \(r_1,r_2,...,r_m\) and \(s_1,s_2,...,s_n\) be sequences of numbers, such that \[ \sum_{k=1}^m r_k=\sum_{k=1}^n s_k \] and \[ \sum_{k=1}^m a_k r_k=\sum_{k=1}^n b_k s_k \] Then \[ \prod_{k=1}^\infty\frac{(k+a_1)^{r_1}(k+a_2)^{r_2}\cdots (k+a_m)^{r_m}}{(k+b_1)^{s_1}(k+b_2)^{s_2}\cdots (k+b_n)^{s_n}}=\frac{\prod_{j=1}^n (b_j!)^{s_j}}{\prod_{j=1}^m (a_j!)^{r_j}} \] In cases where the coefficients are non-integral, we use the Gamma function (an extension of the factorial to non-integers), instead of factorials: \[ \prod_{k=1}^\infty\frac{(k+a_1)^{r_1}(k+a_2)^{r_2}\cdots (k+a_m)^{r_m}}{(k+b_1)^{s_1}(k+b_2)^{s_2}\cdots (k+b_n)^{s_n}}=\frac{\prod_{j=1}^n (\Gamma (b_j+1))^{s_j}}{\prod_{j=1}^m (\Gamma (a_j+1))^{r_j}} \] For instance \[ \prod_{k=0}^\infty \frac{(k+1)(k+a+b)}{(k+a)(k+b)}=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}=B(a,b) \] \[ \frac{\sin(\pi x)}{\pi x}=\prod_{k=1}^\infty\frac{(k-x)(k+x)}{k^2}=\frac{\Gamma(1)^2}{\Gamma(1-x) \Gamma(1+x)}=\frac{1}{\Gamma(1-x)x \Gamma(x)} \] \[ \prod_{k=1}^\infty\frac{(1+\tfrac{1}{k})^x}{1+\tfrac{x}{k}}=\prod_{k=1}^\infty\frac{(k+1)^x k}{k^x (k+x)}=\frac{\Gamma(1)^x \Gamma (1+x)}{\Gamma(2)^x \Gamma(1)}=\Gamma(1+x) \]



Corollary: Asymptotic Behavior of Bernoulli Numbers


In this article, we found that \[ \zeta(2n)=\frac{1}{2}\frac{(2\pi)^{2n}}{(2n)!}\left | B_{2n} \right | \] Combining this with Stirling's approximation, we find that \[ \left | B_{2n} \right |=2\zeta(2n)\frac{(2n)!}{(2\pi)^{2n}} \approx 4\left ( \frac{n}{\pi e} \right )^{2n} \sqrt{n\pi} \cdot e^{1/24n} \]



Corollary: Approximation for Binomial Coefficients


\[ \binom{a}{b}=\frac{a!}{b!(a-b)!} \approx \frac{\left (\tfrac{a}{e} \right )^a \sqrt{2\pi a}} {\left (\tfrac{b}{e} \right )^b \sqrt{2\pi b}\left (\tfrac{a-b}{e} \right )^{a-b} \sqrt{2\pi (a-b)}}=\frac{1}{\sqrt{2\pi}}\sqrt{\frac{a}{b(a-b)}} \frac{a^a}{b^b (a-b)^{a-b}} \]



Corollary: Normal from Binomial


Let \(0 < p < 1\) and \(p+q=1\). Let \[ F_n(x)=\binom{n}{x}p^x q^{n-x} \] \[ f_n(x)=\sqrt{npq}F_n(np+x\sqrt{npq}) \] \[ \phi_n(x)=\ln(f_n(x)) \\ \\ \phi_n(x)=\ln(n!)-\ln((np+x\sqrt{npq})!)-\ln((nq-x\sqrt{npq})!)+(np+x\sqrt{npq})\ln(p)+(nq-x\sqrt{npq})\ln(q) \] Using Stirling's Approximation and some algebra \[ \phi_n(x) = -\tfrac{1}{2}\ln(2\pi)-\left (\tfrac{1}{2}+ np+x\sqrt{npq} \right )\ln\left ( 1+x\sqrt{\frac{q}{np}} \right)-\left (\tfrac{1}{2}+ nq-x\sqrt{npq} \right )\ln\left ( 1-x\sqrt{\frac{p}{nq}} \right)+O(\tfrac{1}{n}) \] Using the series expansion \(\ln(1+x)=x-\tfrac{1}{2}x^2+O(x^3) \) \[ \phi_n(x) = -\tfrac{1}{2}\ln(2\pi)-\tfrac{1}{2}x^2+O(\tfrac{1}{\sqrt{n}}) \] Thus, as \(n\) goes to infinity \[ \phi_\infty(x) = -\tfrac{1}{2}\ln(2\pi)-\tfrac{1}{2}x^2 \] \[ f_\infty(x) = \frac{e^{-x^2/2}}{\sqrt{2\pi}} \] Thus, in the limit, scaling for the changing means and variances, the binomial distribution tends to the normal distribution. Moreover, since the binomial distribution is normalized, we find that \[ \int_{-\infty}^{\infty}\frac{e^{-x^2/2}}{\sqrt{2\pi}}dx=1 \]

Tuesday, October 27, 2015

Occam+Bayes=Induction

  A classic problem in philosophy and the philosophy of science is how to justify induction. That is, how to rationally go from the fact that X is true in N previously observed cases to the belief that it is true in all cases, or at least in an additional, unobserved case. We will here propose a quick and simple method to justify induction, based on the combination of Occam's razor (to choose hypotheses) and Bayesian inference to update epistemic probabilities.


Notation

Let us introduce the following notation. Let \(H\) be some hypothesis which we want to judge for plausibility. Let \(X_k\) be the fact that \(X\) is true in the kth instance. Let \(X^n\) be the fact that \(X\) is true in the first n cases, that is \[X^n=X_1 \cap X_2 \cap \cdots \cap X_n=\bigcap_{k=1}^{n}X_k\] so that \[X^{n-1}\cap X_n=X^n\] Thus \(P\left ( X^n|H \right ) \) is the (epistemic) probability that we observe X in n cases, supposing H is true, and \(P\left ( H|X^n \right ) \) is the (epistemic) probability that H is true, supposing we observe X to be the case in n cases.


Occam's Razor

There are three basic, simplest hypotheses we can form, all the rest being more complex. These three are the
  • Proinductive (P) hypothesis: the chance of X happening again increases as we see more instances of it.
  • Contrainductive (C) hypothesis: the chance of X happening again decreases as we see more instances of it.
  • Uninductive (U) hypothesis: the chance of X happening again stays the same as we see more instances of it.
For concreteness, let \(F_H(n)=P\left ( X_{n}|H \cap X^{n-1} \right )\). Thus we say that, for \(m > 0\), \(F_P(n+m) > F_P(n)\), and \(\lim_{n \rightarrow \infty} F_P(n)=1\), and \(F_C(n+m) < F_C(n)\), and \(\lim_{n \rightarrow \infty} F_C(n)=0\), and \(F_U(n)=F_U(0)\).


Bayesian Inference

We want to find \(P\left ( H|X^n \right ) \) for the hypotheses listed in the previous section. We have \[ P\left ( X^n|H \right )=P\left ( X_n \cap X^{n-1}|H \right )=P\left ( X_n |X^{n-1} \cap H \right ) \cdot P\left ( X^{n-1} |H \right )=F_H(n) \cdot P\left ( X^{n-1} |H \right ) \] Therefore \[ P\left ( X^n|H \right )=\prod_{k=1}^{n} F_H(k) \] Suppose that there are \(N\) mutually exclusive and collectively exhaustive hypotheses. Then, Bayes' formula states: \[ P(H_m|A)=\frac{P(A|H_m)P(H_m)}{P(A|H_1)P(H_1)+P(A|H_2)P(H_2)+\cdots+P(A|H_N)P(H_N)} \] Thus, we have \[ P(H_m|X^n)=\frac{P(X^n|H_m)P(H_m)}{P(X^n|H_1)P(H_1)+P(X^n|H_2)P(H_2)+\cdots+P(X^n|H_N)P(H_N)} \] Therefore \[ P(H_m|X^n)=\frac{P(H_m)\prod_{k=1}^{n} F_{H_m}(k)}{P(H_1)\prod_{k=1}^{n} F_{H_1}(k)+P(H_2)\prod_{k=1}^{n} F_{H_2}(k) + \cdots + P(H_N)\prod_{k=1}^{n} F_{H_N}(k)} \] Let us suppose that the three hypotheses mentioned above are collectively exhaustive. Suppose, for concreteness that \(F_P(n)=\frac{n}{n+1}\), \(F_C(n)=\frac{1}{n+1}\), and \(F_U(n)=\frac{1}{2}\). Thus \(\prod_{k=1}^{n} F_{P}(k)=\frac{1}{n+1}\), and \(\prod_{k=1}^{n} F_{C}(k)=\frac{1}{(n+1)!}\), and \(\prod_{k=1}^{n} F_{U}(k)=\frac{1}{2^n}\). Let \(P(P)=p\) and \(P(C)=q\) and \(P(U)=r\) where \(p+q+r=1\). Then: \[ P(P|X^n)=\frac{p\frac{1}{n+1}}{p\frac{1}{n+1}+q\frac{1}{(n+1)!}+r\frac{1}{2^n}} \] \[ P(C|X^n)=\frac{q\frac{1}{(n+1)!}}{p\frac{1}{n+1}+q\frac{1}{(n+1)!}+r\frac{1}{2^n}} \] \[ P(U|X^n)=\frac{r\frac{1}{2^n}}{p\frac{1}{n+1}+q\frac{1}{(n+1)!}+r\frac{1}{2^n}} \] A simple assessment of limits shows that the former goes to 1 quite rapidly, for increasing n, for any nonzero p, and the latter two go to zero. In fact, for \(p=q=r=1/3\), for \(n>10\), \(P(P|X^n)>0.99\), and for \(n>17\), \(P(P|X^n)>0.9999\).

This example is meant to be only illustrative, to show the general way in which Occam's razor, combined with Bayesian inference, leads to a support of induction. The same things happening repeatedly lends credence to the hypothesis that the same things happen repeatedly, and detracts from the hypothesis that the same things are unlikely to happen repeatedly, or always happen with the same probability. In a very similar way, a coin repeatedly coming up heads supports the hypothesis that it is biased to come up heads, and detracts from the hypotheses that it is biased to come up tails or is fair. This may seem obvious, but it is beneficial to see exactly how the mathematical machinery supports this intuition.

We may also wish to include other hypotheses, but we must first assess the prior probabilities that they are true, and Occam's razor advises taking the inverse probability as inverse to the complexity of the hypothesis. Thus, even if on the hypothesis, observing n X's is more likely than on the three discussed, it would needs be more complex or ad hoc, and so would have a significantly lower prior probability.