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,