InverseErf

*To*: mathgroup at yoda.physics.unc.edu*Subject*: InverseErf*From*: nachbar at merck.com*Date*: Fri, 14 May 1993 11:23:26 -0500 (EDT)

i'm trying to write a package for probability plots, and therefore need to use inverse cumulative distribution functions. the most frequent case involves the normal distribution, which is written in terms or the error function (Erf) in Mathematica. therefore, the inverse functions that i need are written in terms of InverseErf. however, InverseErf seems to be very slow in executing. below is an example (c.f. Out[13]). compiling did not help, see Out[14]. i also tried numerical solutions, which were about 4 times faster (c.f. Out[17])! so what's the bottleneck? requesting Information for the compiled version of my function showed that InverseErf is not compiled. if there is a compiled version in the System` context, how does one get to it? if not, how can one force it to be compiled? could a future release of Mathematica have one? typical applications of probability plot are of the size used in the examples below, and waiting 45 seconds on a unix workstation is too long (i hesitate to try this on my mac iicx!). bob Mathematica 2.1 for SGI Iris Copyright 1988-92 Wolfram Research, Inc. -- Motif graphics initialized -- In[1]:= Needs["Statistics`NormalDistribution`"] In[2]:= Needs["Statistics`InverseStatisticalFunctions`"] In[3]:= Needs["Statistics`DescriptiveStatistics`"] In[4]:= Unprotect[InverseFunction] ; In[5]:= InverseFunction[Erf]=InverseErf ; In[6]:= InverseFunction[Erf,2,2]=InverseErf ; In[7]:= Protect[InverseFunction] ; In[8]:= Solve[p==CDF[NormalDistribution[mu,sigma],x],x] Solve::ifun: Warning: Inverse functions are being used by Solve, so some solutions may not be found. Out[8]= {{x -> mu + Sqrt[2] sigma InverseErf[0, -1 + 2 p]}} In[9]:= InverseCDF[p_,mu_,sigma_]:=N[mu + Sqrt[2] sigma InverseErf[0,2p-1]] In[10]:= CompiledInverseCDF=Compile[{p,mu,sigma},mu + Sqrt[2] sigma InverseErf[0,2p-1]] ; In[11]:= n=100 ; x=Table[(i-.5)/n,{i,n}] ; In[12]:= {m,s}={Mean[x],StandardDeviation[x]} // N Out[12]= {0.5, 0.290115} In[13]:= Timing[InverseCDF[#,m,s]&/@x] Out[13]= {43.75 Second, {-0.247287, -0.129576, -0.0686148, -0.0256623, > 0.00813983, 0.0363403, 0.0607365, 0.0823704, 0.101903, 0.119781,