R/erfinv.R

In LimpopoLab/hydrostats: Statistics for Hydrology

#' Inverse Error Function
#' 
#' This function returns the inverse error function.  The 
#' inverse complementary error function can be obtained 
#' with the arguement of (1-z).  The method is from 
#' Strecok, A. (1968). On the calculation of the inverse 
#' of the error function. Mathematics of Computation<, 
#' 22(101), 144–158. https://doi.org/10.1090/S0025-5718-1968-0223070-2
#' equation 8 with table 2.
#' 
#' @param z argument
#' @return inverse error function
#' @export
erfinv <- function(z) {
      c <- c(0.8862269254527580136490836,
             0.2320136665346544935535340,
             0.1275561753055979582539997,
             0.0865521292415475337296417,
             0.0649596177453854133820146,
             0.0517312819846163741126318,
             0.0428367206517973498446514,
             0.0364659293085316263255797,
             0.0316890050216054468096094,
             0.0279806329649952247334306, #10
             0.0250222758411983494571692,
             0.0226298633188975744328165,
             0.0206067803790590017187687,
             0.0189182172507788544634987,
             0.0174763705628565461904295,
             0.0162315009876852512752949,
             0.0151463150632478055203853,
             0.0141923160025099641511535,
             0.0133473641974212971503471,
             0.0125940048713320698415966, #20
             0.0119182959363920398738582,
             0.0113089701059225367722250,
             0.0107568253033179575777905,
             0.0102542740818534682140955,
             0.0097950057700711756518125,
             0.0093737298191820815168329,
             0.0089859785028433780880376,
             0.0086279535807094346525409,
             0.0082964059277392414044119,
             0.0079885401626033548377731, #30
             0.0077019384322597611663687,
             0.0074344990178315306038419,
             0.0071843865112683124441163,
             0.0069499911010647148571204,
             0.0067298950853415298182763,
             0.0065228451614500617643492,
             0.0063277293643431432300182,
             0.0061435577703841559168944,
             0.0059694462697344562943731,
             0.0058046028538343017603007, # 40
             0.0056483159755515600175138,
             0.0054999446262039635084038,
             0.0053589098416864454702891,
             0.0052246874036874906098435,
             0.0050968015447062055663718,
             0.0049748194997390403988281,
             0.0048583467749584855368178,
             0.0047470230258849766230607,
             0.0046405184555590319253643,
             0.0045385306579070848561814, #50
             0.0044407818435271835322720,
             0.0043470163950215103052998,
             0.0042569987071827288276970,
             0.0041705112741261715711220,
             0.0040873529911090216384896, #55
             0.0040073376434981420567945,
             0.0039302925593064781484431,
             0.0038560574050482175518009,
             0.0037844831074738258156007,
             0.0037154308861260479243302, #60
             0.0036487713836790928372093,
             0.0035843838827445689482347,
             0.0035221555992978687826452,
             0.0034619810441376602308552,
             0.0034037614448720708872919, #65
             0.0033474042218555647207579,
             0.0032928225123033171705625,
             0.0032399347375042310734137,
             0.0031886642086556295315348,
             0.0031389387673655756824422, #70
             0.0030906904573240420296607,
             0.0030438552240413583720698,
             0.0029983726398996133515649,
             0.0029541856520668495355223,
             0.0029112403510908280692453, #75
             0.0028694857582238364863406,
             0.0028288736297367369865527,
             0.0027893582766628551201777,
             0.0027508963985735256089963,
             0.0027134469301298396609447, #80
             0.0026769708992817007406481,
             0.0026414312960977003569962,
             0.0026067929513093085094100,
             0.0025730224237419402268714,
             0.0025400878958849172853698, #85
             0.0025079590769233298461222,
             0.0024766071126183073852077,
             0.0024460045014790969214999,
             0.0024161250167213755471575,
             0.0023869436335520560117534, #90
             0.0023584364613620596151350,
             0.0023305806804456403321513,
             0.0023033544828982969809601,
             0.0022767370173755020142631,
             0.0022507083374217552683530, #95
             0.0022252493531041472147496,
             0.0022003417857069618517663,
             0.0021759681252641061559824,
             0.0021521115907245355018068,
             0.0021287560925625431986507) #100
      su <- 0
      for (n in 1:100) {
            su <- su + c[n] * z^((2*n)-1)
      }
      # test <- abs(erf(su) - z)
      erfinv <- su
      return(erfinv)
}