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$
    – Domen
    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$
    – user1936752
    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$
    – Domen
    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$
    – Domen
    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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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