0
$\begingroup$

Let $\Phi: x\rightarrow y$ be the CDF of a standard normal distribution. The range of $\Phi(x)$ is $[0,1]$ and the domain is all real numbers.

I want to get a series expansion of $\Phi^{-1}(y)$ around $y=0$. I know that the first term should be $-\sqrt{2\ln(y)}$ and would like to know what the higher order terms are.

I tried

f[y_] := InverseCDF[NormalDistribution[0, 1], y] Series[f[y], {y, 0, 5}]

But the result is not a series expansion as expected. There are constant terms like $\log(\frac{1}{2\pi})$ and I'm not sure why they're there since $y\rightarrow 0$. How do I obtain the corrections to the $-\sqrt{2\ln(y)}$ term?

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Use FullSimplify on the result to obtain -Sqrt[-2 Log[y] - Log[-2 π Log[2 π y^2]]] + ... $\endgroup$ Commented Mar 7, 2023 at 21:40
  • $\begingroup$ @Domen thanks! It's a bit annoying that the $\log\log$ term is stuck under the square root - is there a way to separate it from the $\sqrt{-2\log y}$ term? $\endgroup$ Commented Mar 8, 2023 at 12:35
  • $\begingroup$ Well, you can do some funky expression manipulation to remove the other $\log$, like: FullSimplify[Series[f[y], {y, 0, 1}]] /. Log[_ Log[_]] :> 0. But since you mention that you also want to get higher-order terms: I am afraid this is not directly possible in Mathematica. You can see that the expansion is the same, no matter how many terms you specify in Series. That is because Mathematica uses its tabulated series for InverseErfc at $x\to0$, which is a fixed expression. $\endgroup$ Commented Mar 8, 2023 at 22:23
  • 1
    $\begingroup$ Please, take a look at the similar question: How to get more terms with the Series[] expansion of InverseErf[x] around x=1?, and the corresponding answer by @user293787. $\endgroup$ Commented Mar 8, 2023 at 22:24

1 Answer 1

Reset to default
0
$\begingroup$

Ok, mea culpa, I got it wrong.

Try 2:

n = 9; domain = Chop[N[Range[2^(-n), 0.9999999999, 2^(-n)]]]; data = Transpose[{domain, InverseCDF[NormalDistribution[0, 1], domain]}]; indices = (ToString[PaddedForm[#1, 2, NumberPadding -> {"0", "0"}]] & ) /@ Range[0, n]; symbols = (Symbol[StringJoin["x", #1]] & ) /@ indices; expr = symbols . (ChebyshevT[#1, x] & ) /@ Range[0, Length[indices] - 1]; result = FindFit[data, expr, symbols, x]; function = Function[x, Evaluate[FullSimplify[expr /. result]]] 
Function[x, -2.694070753447704 + x*(36.02463367942255 + x*(-414.4544969663766 + x*(2919.8339303927346 + x*(-12190.359149188691 + x*(31174.896173249082 + x*(-49278.882209269956 + x*(46908.437339362164 + x*(-24621.56753111048 + 5471.459451358962*x))))))))] 
Plot[function[y] - InverseCDF[NormalDistribution[0, 1], y], {y, 0, 1}, PlotRange -> {{0, 1}, {-12^(-1), 1/12}}] 
Function[x, -3.1981139183044434 + x*(141.82176399230957 + x*(-7424.74800491333 + x*(236253.08526611328 + x*(-4.680630958251953*^6 + x*(6.148906216479492*^7 + x*(-5.613556024174805*^8 + x*(3.674661709173828*^9 + x*(-1.7557128799768555*^10 + x*(6.145958481669141*^10 + x*(-1.5509398689283594*^11 + x*(2.657566083774026*^11 + x*(-2.4466979730388477*^11 + x*(-9.471151854829346*^10 + x*(6.752732365297246*^11 + x*(-9.872581110395927*^11 + x*(5.544094557729795*^11 + x*(3.4958797462575586*^11 + x*(-8.717278427261528*^11 + x*(6.23268006496337*^11 + x*(-8.207324068464363*^10 + x*(-1.7954943836313046*^11 + x*(1.3475389531949298*^11 + x*(-3.772789388590992*^10 + x*(1.9975620443267794*^9 + 6.922907598366096*^8*x))))))))))))))))))))))))]

The maximum deviation generally stays below $0.05$ until close to the exterme values. The definition of extreme values changes by changes of $n$ values. Try values of $n$ until a satisfactory result is found.

Of course, I still may have missed something.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.