# Asymptotic expansion around infinity for inverse cdf of normal distribution

I'm trying to get a asymptotic expansion as $$x\rightarrow\infty$$ for a particular expression. I have

f[x_] := 1/x*InverseCDF[NormalDistribution[0, 1], 1 - Exp[-x^(1/4)]];


As $$x\rightarrow\infty$$, we have $$f(x)\rightarrow 0$$ but I am interested in how it behaves at $$\infty$$. I tried using AsymptoticSolve and Series but neither seems to get me a nice expansion.

How can I proceed with this? I realize some mathematical insight might also be required before Mathematica can solve this for me.

• Just using Series, Assuming that $x>0$, and doing FullSimplify gives me $\frac{\sqrt{2 \sqrt[4]{x}-\log \left(\pi \left(4 \sqrt[4]{x}-2 \log (2 \pi )\right)\right)}}{x}$. Using series again, one can get the leading term $\sqrt{2} \left(\frac{1}{x}\right)^{7/8}$. More accurately $\sqrt{2} \left(\frac{1}{x}\right)^{7/8}-\frac{\left(\frac{1}{x}\right)^{9/8} \log \left(256 \pi ^4 x\right)}{8 \sqrt{2}}$. Mar 22, 2021 at 8:15
• @yarchik thanks for your comment. Somehow, I still have a giant mess in a square root (and it's different from what's under your square root), even after FullSimplify and Series again. Could you post your code as an answer? Mar 22, 2021 at 8:24

Using assumptions and an immediate assignment:

f[x_] = Assuming[x >= 0,
1/x*InverseCDF[NormalDistribution[0, 1], 1 - Exp[-x^(1/4)]] //
FullSimplify]

(*    -((Sqrt[2] InverseErfc[2 - 2 E^-x^(1/4)])/x)    *)

Assuming[x >= 0, Asymptotic[f[x], x -> ∞]]

(*    Sqrt[2]/x^(7/8) - (8 Log[2] + 4 Log[π] + Log[x])/(
8 Sqrt[2] x^(9/8)) - (
64 Log[2]^2 + 64 Log[2] Log[π] + 16 Log[π]^2 -
64 Log[2 π] + 16 Log[2] Log[x] + 8 Log[π] Log[x] +
Log[x]^2)/(256 Sqrt[2] x^(11/8)) - (1/(
4096 Sqrt[2] x^(
13/8)))(512 Log[2]^3 + 768 Log[2]^2 Log[π] +
384 Log[2] Log[π]^2 + 64 Log[π]^3 -
512 Log[2] Log[2 π] - 256 Log[π] Log[2 π] -
256 Log[2 π]^2 + 192 Log[2]^2 Log[x] +
192 Log[2] Log[π] Log[x] + 48 Log[π]^2 Log[x] -
64 Log[2 π] Log[x] + 24 Log[2] Log[x]^2 +
12 Log[π] Log[x]^2 + Log[x]^3)    *)

Clear["Global*"]

f[x_] = Assuming[x >= 0,
1/x*InverseCDF[NormalDistribution[0, 1], 1 - Exp[-x^(1/4)]] // FullSimplify]

(* -((Sqrt[2] InverseErfc[2 - 2 E^-x^(1/4)])/x) *)

f1[x] = Asymptotic[f[x], x -> ∞] // Simplify[#, x >= 0] &

(* Sqrt[2 x^(1/4) - Log[2 π (2 x^(1/4) + Log[1/(2 π)])]]/x *)

f2[x] = Assuming[x >= 0, Asymptotic[f[x], x -> ∞] // FullSimplify]

(* (64 (128 x^(3/4) + 4 Log[2 π]^2 - 32 Sqrt[x] Log[4 π] +
4 Log[2 π] Log[4 π] - Log[4 π]^3 -
4 x^(1/4) (4 Log[2]^2 - 4 Log[2 π] + Log[π] Log[16 π])) -
Log[x] (16 (32 Sqrt[x] - 4 Log[2 π] + 8 x^(1/4) Log[4 π] +
3 Log[4 π]^2) + 4 (4 x^(1/4) + Log[64 π^3]) Log[x] +
Log[x]^2))/(4096 Sqrt[2] x^(13/8)) *)


Plotting the functions, f2 appears to be a better approximation

Plot[Evaluate@{f1[x], f2[x], f[x]}, {x, 5, 40}, PlotLegends ->
Placed[{"f1", "f2", "f"}, {.5, .5}]]


However, looking at the relative error, f1 is better for larger values of x

LogPlot[Evaluate@{(f1[x] - f[x])/f[x], (f2[x] - f[x])/f[x]}, {x, 5, 150},
WorkingPrecision -> 50,
PlotLegends ->
Placed[{"f1", "f2"}, {.75, .75}]]


The crossover occurs at

FindRoot[f1[x] == f2[x], {x, 90}]

(* {x -> 91.366} *)


a Piecewise` function could be used to switch between the approximations.