I programmed up a little algorithm that goes like this:

  1. Fix two positive, real numbers, call them $\alpha$ and $\beta$.
  2. Generate a new, random, real number, $x \in [0,1]$
  3. Set $\alpha$ = $x\alpha\beta$,
  4. Repeat 2 and 3 until $\alpha$ is equal to the smallest float supported (I'll just call this 0), or $\alpha$ is equal to the largest float supported (I'll call this $\infty$).

Initially I thought this would turn up nonsense, or only occasionally go to 0 or infinity.

But after experimenting a lot with values of $\beta$, I found a fascinating result:

If $\beta \in (0,e)$, $\alpha$ always goes to $0$.

If $\beta=e$, $\alpha$ goes to $0$ about 50% of the time, and to $\infty$ about 50% of the time.

If $\beta \in (e,\infty)$, $\alpha$ always goes to $\infty$.

Now these are all results based on running the above algorithm about 10,000 times, so maybe I'm in the wrong to put 'always', but certainly most of the time these results hold.

This was surprising to me! I'm not well versed in randomness, or computational mathematics, or probability so I can't really fathom why $e$ would pop up so prominently in this context. Can anyone point me towards the subject I need to examine to understand this result?

(Apologies if this post is mistagged, as I said, I'm not experienced with the subjects this post related to)

  • $\begingroup$ The implementation I've written is in Ruby, using the Random.rand method. I'm not familiar with PRNGs, but documentation says "PRNGs are currently implemented as a modified Mersenne Twister with a period of 2**19937-1." This all generates a random float between 0 and 1, which I then multiply by $\beta$. $\endgroup$ – Emma Foster Mar 10 '14 at 5:15
  • 1
    $\begingroup$ I think you are multiplying by a number between $0$ and $\beta$ at each stage, right? If so, you need to either write "$x\in[0,1]$" in step 2 or remove the factor of $\beta$ in step 3. $\endgroup$ – mjqxxxx Mar 10 '14 at 6:51
  • $\begingroup$ Yes, you're correct. I've modified step 2. $\endgroup$ – Emma Foster Mar 10 '14 at 7:06
  • 1
    $\begingroup$ The crucial insight is that if a random variable $X$ is chosen from a uniform distribution from $[0, 1]$, then $-\ln X$ is from an exponential distribution with mean $1$. (And this is how $e$ pops up.) And so, when you multiply by $\beta$, it matters whether $\ln \beta$ has mean more than $1$ or not; this determines whether it offsets $\ln X$ with mean $-1$ or not. $\endgroup$ – ShreevatsaR Mar 11 '14 at 1:29

For a fixed $\beta$, your algorithm can be stated as follows:

  • Pick $\alpha_0$ arbitrarily.
  • For each $n \ge 1$, generate a uniformly random real number $X_n$ in $[0, 1]$ and set $\alpha_n = X_n\alpha_{n-1}\beta$.

Thus each $\alpha_n$ is got by multiplying the previous $\alpha_{n-1}$ by $\beta X_i$, so $$\alpha_n = \alpha_0\beta^n X_1 X_2 \cdots X_n \tag 1$$

Things become more tractable/familiar if we take logarithms, to convert products into sums: each $\log \alpha_n$ is got by taking the previous $\log \alpha_{n-1}$ and adding $\log\beta + \log X_i$, so we get the logarithmic version of $(1)$: $$\log \alpha_n = \log \alpha_0 + n\log\beta + \sum_{i=1}^n \log X_i \tag 2$$

Thus the sequence $\log \alpha_n$ is a random walk starting at $\log \alpha_0$ and changing by $Y_i = \log \beta + \log X_i$ at each step. For such a random walk $\sum_{i=1}^n Y_i$, its eventual behaviour is determined by whether $E[Y_i]$ is positive, negative, or zero.

We can in fact determine the distribution of $\log X_i$ (and therefore of $Y_i$). Each $X_i$ is uniformly distributed in $[0, 1]$. It turns out, this means that $(-\ln X_i)$ is exponentially distributed with mean $1$. Let's derive that. To find the distribution of $\log_b X_i$ (taking logarithms to arbitrary base $b$ to be clear on where $e$ pops up), we we can calculate the cumulative distribution function: for $t \le 0$ (assuming wlog that $b > 1$), we have $\Pr(\log_b X_i \le t) = \Pr(X_i \le b^t) = b^t$. Differentiating this wrt $t$ gives the probability density function, which is $b^t\ln b$ (and $0$ for $t > 0$). In particular, $$ \operatorname{E}[\log_b X_i] = \int_{-\infty}^{0} t \, b^t\ln b \, dt = - 1/\ln b = -\log_b e. $$

This is where $e$ pops up: although to convert products into sums we could have taken logarithm to any base $b$, the fact is that $E[\log_b X_i] = -\log_b e$, so $e$ pops up unavoidably (and we might as well take the natural logarithm).

So $E[Y_i] = \log\beta + \log X_i = \log \beta - \log e$ (whatever base we choose), which is positive, negative or zero depending on whether $\beta > e$, $\beta < e$ or $\beta = e$. It is therefore intuitively clear that when $\beta > e$ the random walk $\log\alpha_n$ drifts towards $\infty$ (and $\alpha_n$ towards $\infty$ as well) and when $\beta < e$ towards $-\infty$ (and $\alpha_n$ towards $0$). Proofs of precise statements below.

Let $S_n = \sum_{k=1}^n \log X_i$. We have, from the strong law of large numbers applied to the independent and identically distributed random variables $\log X_i$, the fact ("a.s." stands for "almost surely") that $$\frac{S_n}{n} \xrightarrow{a.s.} - 1.$$

This means that, with probability $1$, we have $$\frac{\log \alpha_n}{n} = \frac{\log \alpha_0}{n} + \log \beta + \frac{S_n}{n} \xrightarrow{a.s.} \log \beta - 1$$

In your problem, we iterate until either $\alpha_n$ hits an upper threshold or hits a lower threshold, and we record which threshold it hit. (Typical example thresholds induced by programming languages are on the order of $10^{38}$ and $10^{-38}$ for single-precision floating-point format, and $10^{308}$ and $10^{-308}$ for double-precision floating-point format.) In particular,

  • when $\beta < e$, we have $\lim\limits_{n\to\infty}\dfrac{\log \alpha_n}{n} < 0$ almost surely, which means that $a_n \to 0$ almost surely (and thus it is more likely to hit the upper threshold),

  • when $\beta > e$, we have $\lim\limits_{n\to\infty}\dfrac{\log \alpha_n}{n} = c > 0$ almost surely (where $c = \log\beta - 1$), which means that $a_n \to \infty$ almost surely ($\alpha_n$ grows like $e^{nc} \to \infty$) (and thus it is more likely to hit the upper threshold).

For the case when $\beta = e$, see below.

Note however that as our thresholds are not exactly $0$ and $\infty$ (neither of which is achievable anyway), we can only say "more likely" above: the fact that almost surely $a_n \to 0$ doesn't necessarily mean that the sequence $a_n$ hits $10^{-308}$ before $10^{308}$: it is possible (just very unlikely) that it does reach $10^{308}$ before $10^{-308}$, but eventually approaches $0$.

I don't know how to fix this, but at least for the $\beta = e$ case, the theory of random walks (see section 9.4) has some help. Let $Y_i = \log\beta + X_i$, and let $T_n = \sum_{i=0}^{n} Y_i$ (using $T_n$ to avoid confusion with the earlier $S_n$) so that $T_n$ is a random walk, and $\log\alpha_n = \log\alpha_0 + T_n$. Then if the lower threshold is $L$ and the upper threshold is $U$, the stopping events $\alpha_n \le L$ and $\alpha_n \ge U$ are equivalent to $T_n \le \log L - \log \alpha_0 = L'$ and $T_n \ge \log U + \log\alpha_0 = U'$.

For the special case when $\operatorname{E}[Y_i] = 0$ (this means $\beta = e$), we have, denoting $J$ as the stopping time, $\operatorname{E}[S_J] = \operatorname{E}[J]\operatorname{E}[Y_i] = 0$ (Wald's equality), which roughly gives, assuming that $S_J$ is either roughly $L'$ or $U'$ (i.e., $\alpha_J$ doesn't overshoot the thresholds by too much) that $0 = \operatorname{E}[S_J] = p_{L'}L' + (1-p_{L'})U'$ so $p_{L'} = \frac{U'}{U'-L'} = \frac{\log U + \log\alpha_0}{\log U - \log L + 2\log \alpha_0}$ which for the thresholds from the single-precision floating format $U \approx 3.4 × 10^{38}$ and $L = 2^{−149}$ gives roughly $\frac{88.7 + \log \alpha_0}{192.0 + 2\log\alpha_0}$ which at least for values of $\log\alpha_0$ small in magnitude gives roughly $0.46$ rather than $0.5$ exactly (though it's much closer to $0.50$ if you take your $L$ to be the smallest normalized number $2^{-126}$ instead).

  • 1
    $\begingroup$ @TamTamTam: Actually I think there is more to say here... look back in a day or so; there may be more answers here (or even a bit more in my answer :p). $\endgroup$ – ShreevatsaR Mar 10 '14 at 5:39
  • 2
    $\begingroup$ If $\alpha_{n}=X_{n}\alpha_{n-1}\beta$, then $\alpha_{n}=\alpha_{0}\beta X_1 \beta X_2 \cdots \beta X_n$. You are missing many $\beta$'s. $\endgroup$ – mjqxxxx Mar 10 '14 at 6:19
  • 1
    $\begingroup$ My code that simulates this (with sample output) is here: github.com/shreevatsa/misc-math/blob/master/log-convergence.py $\endgroup$ – ShreevatsaR Mar 10 '14 at 7:17
  • 1
    $\begingroup$ I do think there is still something interesting going on, which is that the simulation doesn't show a sharp transition around $\beta = e$. $\endgroup$ – ShreevatsaR Mar 10 '14 at 7:21
  • 1
    $\begingroup$ I think that it is OK now. Your code's behaviour around $e$ also interests me, so I will give it a try in an online compiler while on cellphone, and will comment here the result if I survive to such an adventure :-) $\endgroup$ – chubakueno Mar 10 '14 at 7:32

You have $\alpha_{i+1}=\alpha_{i} X_{i+1}$, where $X_i\sim U(0,\beta)$; so $\alpha_{n}=\alpha_{0}\prod_{i=1}^{n}X_{i}$. The infinite product (and hence $\alpha_n$) diverges to $+\infty$ if its logarithm, $$ \sum_{i=1}^{n}\log X_i, $$ does, and by the law of large numbers this occurs almost surely if $E[\log X_i]>0$. Now, $$ E[\log X_i]=\frac{1}{\beta}\int_{0}^{\beta}\log x dx=\log(\beta)-1=\log(\beta/e), $$ which is positive when $\beta>e$. Similarly, the product diverges to $0$ if its logarithm diverges to $-\infty$, which occurs almost surely if $\beta < e$. Finally, if $\beta=e$ exactly, you have an unbiased one-dimensional random walk, which implies that $\alpha_n$ neither converges nor diverges: with probability $1$, it crosses $1$ (that is, goes from larger than $1$ to smaller, or vice-versa) infinitely many times.

  • 1
    $\begingroup$ OP's observation… that the cutoff is $\beta=e$… can only be correct if there's something wrong with how the question is written. My guess is that the recursion is actually $\alpha_{i+1}=\alpha_{i}x_{i+1}$; i.e., the extra "$\beta$" is extraneous. $\endgroup$ – mjqxxxx Mar 10 '14 at 5:28
  • $\begingroup$ No, it shouldn't. $X$ is uniformly chosen from $[0,\beta^2]$. $\endgroup$ – mjqxxxx Mar 10 '14 at 6:25
  • $\begingroup$ This looks right +1 - my code shows somewhere $\beta \approx 1.65$ as the cutoff. Will think more on this. $\endgroup$ – Macavity Mar 10 '14 at 6:31
  • $\begingroup$ @Macavity: Was about to leave this comment on your deleted answer: By the way, why do you run it for a fixed number NMax of iterations, instead of until it hits a threshold (as in the question)? And are you visually inspecting the value of a at the end of your NMax iterations, or averaging it over different runs (in which case the expected value will become relevant)? $\endgroup$ – ShreevatsaR Mar 10 '14 at 6:56
  • $\begingroup$ @ShreevatsaR and Macavity can you please share your full source code? I would like to contrast them and see if I find something.(I don't have a particular simulation result to share, as I am using phone right now ) $\endgroup$ – chubakueno Mar 10 '14 at 7:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.