Asymptotics of a randomized Fibonacci sequence Let $f(1)=f(2)=1$ and recursively define $f(n+1) = f(n) + f(i)$, where $i$ is chosen uniformly at random from $1,\ldots,n-1$. About how big should we expect $f(n)$ to be for $n$ large? We can examine the expected value of $f(n)$, say $E(f(n)) = g(n)$. After some algebra we get the second-order difference equation $$g(n+2) = \left(2-\frac{1}{n}\right)g(n+1) - \left(1-\frac{2}{n}\right)g(n)$$ for $n>0$. I did some guesswork by first looking at related differential equations, then adjusting the fit, then numerically estimating the constant to suggest $$g(n) \sim \frac{1}{\sqrt{\pi e}} \frac{e^{2\sqrt{n}}}{n^{3/4}}.$$
That the order is correct follows from Wong and Li, "Asymptotic Expansions for Second-Order Linear Difference Equations".
I have two questions.
First question. How do you rigorously prove that $1/\sqrt{\pi e}$ is the correct proportionality constant? I guessed it from numerical evidence, but the fit is good enough that I believe it's correct. The paper by Wong and Li doesn't appear to address determining the proportionality constant, likely because it depends on initial values.
Second question. Wong and Li give a recursive procedure at (2.18) to generate higher-order asymptotic terms. I've implemented the code to compute the next term $c_1$, but it disagrees with the value for $c_1$ they give in an alternate form at (2.19). Numerically, the value from (2.19) appears to be correct. Is there a typo in (2.18) somewhere?
Using (2.19) we get $$g(n) \sim \frac{1}{\sqrt{\pi e}} \frac{e^{2\sqrt{n}}}{n^{3/4}}\left(1-\frac{101}{48}\frac{1}{\sqrt{n}}\right).$$
I've implemented both (2.18) and (2.19) to compute $c_1$ in SageMath (my code).
My SageMath code calculates $187/48$ using (2.18), but $-101/48$ using (2.19). The value $-101/48$ appears to be correct numerically.
Added: Brendan McKay pointed out this wonderful connection to sequence A000262 in the OEIS: $$g(n) = \frac{2}{(n-2)!} A(n-2)$$ for $n>2$. I haven't proven this relationship, but it definitely appears to be true. It implies, for example, that $$g(n) = \frac{2}{n-2} L_{n-3}^1 (-1),$$ the generalized Laguerre polynomials, for $n>2$, giving the asymptotic series $$g(n) = 
\frac{1}{\sqrt{\pi e}} \frac{e^{2\sqrt{n}}}{n^{3/4}}\left(
1 - \frac{101}{48}\frac{1}{n^{1/2}}
+ \frac{16609}{4608}\frac{1}{n}
- \frac{18575473}{3317760}\frac{1}{n^{3/2}}
+ \ldots
\right).$$
This establishes the proportionality constant for the asymptotics of this particular sequence, but still leaves open my original questions. How do you determine proportionality constants for similar second-order linear difference equations in general? Given the additional terms for this asymptotic series, I should be able to track down the error(s) in the recursive algorithm given by Wong and Li, but it'll take some calculation and time.
Added: I've now proven that $g(n) = 2 L_{n-2}^{-1}(-1)$, and it was actually quite easy.
We have the recurrence $$g(n+3) = \left(2-\frac{1}{n+1}\right) g(n+2) - \left(1-\frac{2}{n+1}\right) g(n+1),$$ which we can write $$g(n+3) = \frac{1}{n+1} \left( (2n+1) g(n+2) - (n-1) g(n+1) \right).$$
But this is a generalized Laguerre recurrence!
We need $1+\alpha-x = 1$ and $-1 = \alpha$, and so $\alpha = x = -1$, thus $$g(n) = C\cdot L_{n-2}^{-1}(-1)$$ for some constant $C$, assuming agreement with initial values. But $g(3) = 2 = 2 L_1^{-1}(-1) = 2\cdot 1$ and $g(4) = 3 = 2 L_2^{-1}(-1) = 2\cdot 3/2$, thus $g(n) = 2 L_{n-2}^{-1}(-1)$ for all $n>2$!
 A: (this needs to be updated for the correction above)
You can compute the generating function via:  Let $M(z) = \sum_{n=1} g(n+1) z^n$
Then the recursion relation for $g(n+1)$ gives the relation
$$M(z) = g(1) + z g(2) +$$$$ +2 z (M(z) - g(1)) - \int ( M(z) - g(1)) dz + 2 \int ( z M(z) ) dz - z^2 M(z)$$
differentiate both sides of this relation to eliminate the anti-derivative terms:
$$M'(z) = g(2) + 2 (M(z) - g(1)) +$$$$+ 2 z M'(z) - M(z) + g(1) + 2 z M(z) - 2 z M(z) - z^2 M'(z)$$
solve for $M'(z)$:
$(z^2 - 2 z + 1) M'(z) = M(z)$ (using the initial condition simplifies things)
and then we have
$$d/dz \log(M(z)) = 1/(1 - z)^2 $$
which is solved by
$$M(z) = K e^{1/(1-z)}$$
Using the initial condition again we find $k = 1/e.$
Fundamental singularity of $z=1$ giving a radius of convergence of $1$ will determine asymptotics of the coefficients. Exact formulas for the $g(n)$ should also be possible.
A: The question has been solved with an appeal to Laguerre polynomials, but for a more general approach, you can derive a generating function like that done by contributor Virgil.  Then asymptotics with an appropriate leading factor can be derived with variations of the saddle point method.  In 'Asymptotics of Decomposable Combinatorial Structures of Alg-Log Type with Positive Log Exponent,' by Z. Gao, et. al., a derivation for the asymptotics of the coefficients of an exponential of a generating function $G(z)$ that is of rather general functional form:
$$ \text{For } G(z) \sim(1-z)^{-a}\Big( \log(1/(1-z)) \Big) ^b ,\quad a>0, b>0$$
$$ \text{Let } c=a^{-1/(1+a)}(1+a)^{1+b\frac{a+2}{a+1}}, r\sim 1-(a/n)^{1/(1+a)}\Big( \frac{\log(n)}{1+a} \Big)^{b/(1+a)} $$
$$ \text{Then  }\quad
[z^n] \exp{(G(z))}  \sim \exp{(G(r))}r^{-n}\sqrt{\frac{1}{2\pi\ c} n^{-\frac{a+2}{a+1}} \big(\log{n} \big)^{\frac{b}{1+a}}   }
$$
This does not apply to your situation because $b=0.$  However, it illustrates the generality of the result.  The derivation uses Hayman's method. The paper also treats $a=0, b \ge 2.$  References are given for other cases as well.
