Julia に翻訳--220 ppnd，正規分布のパーセント点（0.975 を与えると 1.96 を返すような）

#==========
Julia の修行をするときに，いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。

algorithm as241  appl. statist. (1988) vol. 37, no. 3, 477- 484.
https://jblevins.org/mirror/amiller/as181.f90

ファイル名: ppnd.jl　　関数名: ppnd

翻訳するときに書いたメモ

元は FORTRAN プログラム。書き方も古風なので（GO TO があったり），少し新しめに。

以下の処理を明確化
p = 0, 1 ==> -Inf, Inf
p < 0, p > 1 ==> NaN

==========#

function ppnd(p)
  #=====
  algorithm as241  appl. statist. (1988) vol. 37, no. 3, 477- 484.

  produces the normal deviate z corresponding to a given lower tail area of p;
  z is accurate to about 1 part in 10^7.
  =====#

# coefficients for p close to 0.5

  a0 = 3.3871327179; a1 = 50.434271938; a2 = 159.29113202; a3 = 59.10937472
  b1 = 17.895169469; b2 = 78.757757664; b3 = 67.1875636

# coefficients for p not close to 0, 0.5 or 1.

  c0 = 1.4234372777; c1 = 2.7568153900; c2 = 1.3067284816; c3 = 0.17023821103
  d1 = 0.7370016425; d2 = 0.12021132975

# coefficients for p near 0 or 1.

  e0 = 6.6579051150; e1 = 3.0812263860; e2 = 0.42868294337; e3 = 0.017337203997
  f1 = 0.24197894225; f2 = 0.012258202635

  q = p - 0.5
  if abs(q) <= 0.425
    r = 0.180625 - q^2
    z = q * (((a3 * r + a2) * r + a1) * r + a0) / (((b3 * r + b2) * r + b1) * r + 1)
  else
    r = q < 0 ? p : 1 - p
    if r < 0
      z = NaN
    elseif r == 0
      z = Inf
    else
      r = sqrt(-log(r))
      if r <= 5
        r -= 1.6
        z = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + 1)
      else
        r -= 5
        z = (((e3 * r + e2) * r + e1) * r + e0) / ((f2 * r + f1) * r + 1)
      end
    end
    q < 0 && (z = -z)
  end
  z
end

R では qnorm() に相当する。

using Rmath
function check(x)
  println("ppnd($x) = $(ppnd(x)), qnorm($x) = $(qnorm(x)), abs(diff) = $(abs(ppnd(x)-qnorm(x)))")
end

check.([0.975, 0.0001, 0.4, -1, 0, 1, 2])
# ppnd(0.975) = 1.9599641740135167, qnorm(0.975) = 1.9599639845400536, abs(diff) = 1.8947346314135416e-7
# ppnd(0.0001) = -3.719016535339939, qnorm(0.0001) = -3.71901648545568, abs(diff) = 4.988425894580928e-8
# ppnd(0.4) = -0.2533470936835343, qnorm(0.4) = -0.2533471031357998, abs(diff) = 9.452265470333288e-9
# ppnd(-1.0) = NaN, qnorm(-1.0) = NaN, abs(diff) = NaN
# ppnd(0.0) = -Inf, qnorm(0.0) = -Inf, abs(diff) = NaN
# ppnd(1.0) = Inf, qnorm(1.0) = Inf, abs(diff) = NaN
# ppnd(2.0) = NaN, qnorm(2.0) = NaN, abs(diff) = NaN