Julia に翻訳--223 Shapiro-Wilk 検定，正規性の検定

Julia に翻訳--223　Shapiro-Wilk 検定，正規性の検定

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

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

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

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

Julia に翻訳--220, 221, 222 は他ならず，これを翻訳するための道具であった。
alnorm, ppnd は名前をそのまま残し，pnorm, qnorm を使う。
alnorm(z, flag=true) = pnorm(z, flag)
ppnd(p) = qnorm(p)
途中，謎の数値がいくつかあるが，正確な記述に改めた。
結果として R の shapiro.test() と同じ精度の結果が得られる。（n = 3 のときのp値の末尾がちょっと違う）
その他，Julia の特徴を生かした簡潔なプログラムになったと思う。

==========#

using Statistics, Rmath

function shapirotest(x)
  x = sort(x);
  n = length(x)
  (n < 3 || n > 5000) && error("sample size must be between 3 and 5000")
  rng = x[n] - x[1]
  rng == 0 && ("all 'x' values are identical")
  rng < 1e-10 && (x ./= rng)
  # 仮定
  # 1. 打切りデータ，欠損値はない
  # 2. 3 <= n <= 5000
  # 3. データは昇順にソートされている
  # 4. 全て同じ値ということはない
  # 5. rng < 1e-10 のときは，事前に標準化されている
  swilk(x, n)
end

function swilk(x, n)
  #=====
           algorithm as r94 appl. statist. (1995) vol.44, no.4

           calculates the shapiro-wilk w test and its significance level

    input:
      x(n)    sample values in ascending order.
      n        the total sample size (including any right-censored values).
    output
      w        the shapiro-wilks w-statistic.
      pw       the p-value for w.
  =====#

  n2 = n ÷ 2

  # calculates coefficients for the test
  if n == 3
    # a = [0.70711]
    a = [1/sqrt(2)]
  else
    a = zeros(n2);
    an25 = n + 0.25
    a = [ppnd((i - 0.375) / an25) for i = 1:n2]
    summ2 = 2sum(a .^ 2)
    ssumm2 = sqrt(summ2)
    rsn = 1 / sqrt(n)
    c1 = [0.0, 0.221157, -0.147981, -2.07119, 4.434685, -2.706056]
    a1 = poly(c1, rsn) - a[1] / ssumm2

    # normalize coefficients
    if n > 5
      c2 = [0.0, 0.042981, -0.293762, -1.752461, 5.682633, -3.582633]
      a2 = -a[2] / ssumm2 + poly(c2, rsn)
      fac = sqrt((summ2 - 2 * a[1] ^ 2  - 2 * a[2] ^ 2) / (1 - 2 * a1 ^ 2 - 2 * a2 ^ 2))
      a[1] = a1
      a[2] = a2
      a[3:n2] ./= -fac
    else
      fac = sqrt((summ2 - 2 * a[1] ^ 2) / (1 - 2 * a1 ^ 2))
      a[1] = a1
      a[2:n2] ./= -fac
    end
  end

  # calculate w statistic as squared correlation between data and coefficients
  x ./= x[n] - x[1]
  x .-= mean(x)
  ssa = ssx = sax = 0
  for i = 1:n
    j = n + 1 - i
    asa = i == j ? 0 : sign(i - j) * a[min(i, j)]
    ssa += asa^2
    ssx += x[i]^2
    sax += asa * x[i]
  end

  # w1 equals (1-w) claculated to avoid excessive rounding error
  # for w very near 1 (a potential problem in very large samples)
  ssassx = sqrt(ssa * ssx)
  w1 = (ssassx - sax) * (ssassx + sax) / (ssa * ssx)
  w = 1 - w1

  # calculate significance level for w (exact for n=3)
  if n == 3
    # pw = 1.909859 * (asin(sqrt(w)) - 1.047198)
    pw = 6 / π * (asin(sqrt(w)) - π / 3)
  else
    y = log(w1)
    xx = log(n)
    m = 0
    s = 1
    if n <= 11
      gamma = poly([-2.273, 0.459], n)
      y >= gamma && return w, 1e-10
      y = -log(gamma - y)
      m = poly([0.5440, -0.39978, 0.025054, -0.6714e-3], n)
      s = exp(poly([1.3822, -0.77857, 0.062767, -0.0020322], an))
    else
      m = poly([-1.5861, -0.31082, -0.083751, 0.0038915], xx)
      s = exp(poly([-0.4803, -0.082676, 0.0030302], xx))
    end
    pw = alnorm((y - m) / s, false)
  end
  w, pw
end

function poly(coefficients, x)

  #        algorithm as 181.2   appl. statist.  (1982) vol. 31, no. 2

  #        calculates the algebraic polynomial of order nored-1 with
  #        array of coefficients coefficients.  zero order coefficient is coefficients[1]

  nord = length(coefficients)
  result = 0
  for i = nord:-1:2
    result = (result + coefficients[i]) * x
  end
  result + coefficients[1]
end

ppnd(p) = qnorm(p)
alnorm(z, flag=true) = pnorm(z, flag)

使用例

偶数個
x = [ 0.03, 0.08, 0.14, 0.14, 0.15, 0.17, 0.31, 0.37, 0.37, 0.40,
      0.42, 0.43, 0.49, 0.53, 0.60, 0.62, 0.72, 0.75, 0.87, 0.97 ];

奇数個
y = [ 0.03, 0.08, 0.14, 0.14, 0.15, 0.17, 0.31, 0.37, 0.37, 0.40,
      0.42, 0.43, 0.49, 0.53, 0.60, 0.62, 0.72, 0.75, 0.87 ];

最小個数
z = [ 0.03, 0.08, 0.14 ];

最大個数
using Random
Random.seed!(321)
u = randn(5000);

shapirotest(x) # (0.9594825287257593, 0.5335636425303667)
shapirotest(y) # (0.9607205920932019, 0.5866579913499355)
shapirotest(z) # (0.9972527472527474, 0.8998502800372317)
shapirotest(u) # (0.9995093332208211, 0.22746134654936953)

R の shapiro.test() の結果と比較
using RCall
R"""
options(digits=16)
a = shapiro.test($x); cat(a$statistic, a$p.value,"\n")
b = shapiro.test($y); cat(b$statistic, b$p.value,"\n")
c = shapiro.test($z); cat(c$statistic, c$p.value,"\n")
d = shapiro.test($u); cat(d$statistic, d$p.value,"\n")
"""
# 0.9594825287257593 0.5335636425303667
# 0.9607205920932019 0.5866579913499355
# 0.9972527472527474 0.8998502800372251
# 0.9995093332208211 0.2274613465493695

Julia に翻訳--222 多項式の値を求める

Julia に翻訳--222　多項式の値を求める

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

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

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

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

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

==========#

function poly(coefficients, x)

  #        algorithm as 181.2   appl. statist.  (1982) vol. 31, no. 2

  #        calculates the algebraic polynomial of order nored-1 with
  #        array of coefficients coefficients.  zero order coefficient is coefficients[1]

  nord = length(coefficients)
  result = coefficients[1]
  nord == 1 && return result
  p = x * coefficients[nord]
  if nord != 2
    n2 = nord - 2
    j = n2 + 1
    for i = 1:n2
      p = (p + coefficients[j]) * x
      j -= 1
    end
  end
  result + p
end

元のプログラムは古い FORTRAN で書かれているので，もっとわかりにくいが，書き換えてもなお，
プログラムが長く，一見何をやっているかわかりづらいが，nord が 1 であるかとか，2 であるとか，定数項を特別扱いしたりとかしなくてもよい。
たとえば，以下の5次式
0.1 + 0.221157∙x - 0.147981∙x^2 - 2.07119∙x^3 + 4.434685∙x^4 - 2.706056∙x^5
を ^ を使わないで
((((-2.706056∙x + 4.434685)∙x - 2.07119)∙x - 0.147981)∙x + 0.221157)∙x + 0.1
として計算する。

以下のように，短く簡潔でわかりやすいプログラムになる。

function poly2(coefficients, x)

  #        algorithm as 181.2   appl. statist.  (1982) vol. 31, no. 2

  #        calculates the algebraic polynomial of order nored-1 with
  #        array of coefficients coefficients.  zero order coefficient is coefficients[1]

  nord = length(coefficients)
  result = 0
  for i = nord:-1:2
    result = (result + coefficients[i]) * x
  end
  result + coefficients[1]
end

c6 = [0.1, 0.221157, -0.147981, -2.07119, 4.434685, -2.706056];
# 0.1 + 0.221157∙x - 0.147981∙x^2 - 2.07119∙x^3 + 4.434685∙x^4 - 2.706056∙x^5
# ((((-2.706056∙x + 4.434685)∙x - 2.07119)∙x - 0.147981)∙x + 0.221157)∙x + 0.1 これを計算する

x = 1.2
poly(c6, x)  # -0.9644910099199991
poly2(c6, x) # -0.9644910099199991

Julia では以下のようにすればよい。
using Polynomials
func = Polynomial(c6)
typeof(func) # Polynomial{Float64, :x}; func は関数
func(x)      # -0.9644910099199991

2次式，1次式，0次式の場合もちゃんと動くのを確認しておく。

c3 = [1.2, 3.4, 5.6]; # (5.6∙x + 3.4)∙x + 1.2
poly(c3, 3)       # 61.79999999999999
poly2(c3, 3)      # 61.79999999999999
Polynomial(c3)(3) # 61.79999999999999

c2 = [1.2, 3.4]; # 3.4∙x + 1.2
poly(c2, 3)       # 11.399999999999999
poly2(c2, 3)      # 11.399999999999999
Polynomial(c2)(3) # 11.399999999999999

c1 = [1.2]; # 1.2; x に何を指定しても無関係（定数が返される）
poly(c1, 5)       # 1.2
poly2(c1, 5)      # 1.2
Polynomial(c1)(5) # 1.2

Julia に翻訳--221 alnorm，正規分布の確率（1.96 を与えると 0.975 を返すような）

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

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

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

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

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

==========#

function alnorm(z, lower=true)

    #=====
       algorithm as66 applied statistics (1973) vol.22, no.3

       evaluates the tail area of the standardised normal curve
       from x to infinity if upperper is true or
       from minus infinity to x if upperper is .false.
    =====#

    P = 0.398942280444; Q = 0.39990348504; R = 0.398942280385;
    a1 = 5.75885480458; a2 = 2.62433121679; a3 = 5.92885724438
    b1 = -29.8213557807; b2 = 48.6959930692
    c1 = -3.8052e-8; c2 = 3.98064794e-4; c3 = -0.151679116635
    c4 = 4.8385912808;c5 = 0.742380924027; c6 = 3.99019417011
    d1 = 1.00000615302; d2 = 1.98615381364; d3 = 5.29330324926
    d4 = -15.1508972451; d5 = 30.789933034

    z > 0 && (lower = !lower)
    z = abs(z)
    if z <= 7 || !lower && z <= 18.66
        y = 0.5 * z^2
        if z > 1.28
            p = R * exp(-y) / (z + c1 + d1 / (z + c2 + d2 /
                (z + c3 + d3 / (z + c4 + d4 / (z + c5 + d5 / (z + c6))))))
        else
            p = 0.5 - z * (P - Q * y / (y + a1 + b1 / (y + a2 + b2 / (y + a3))))
        end
    else
        p = 0
    end
    !lower && (p = 1 - p)
    p
end

alnorm(1.96)         # 0.9750021048508247
alnorm(1.96, false)  # 0.024997895149175307
alnorm(-1.96)        # 0.024997895149175307
alnorm(-1.96, false) # 0.9750021048508247

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

using Rmath
pnorm(1.96)          # 0.9750021048508247
pnorm(1.96, false)   # 0.024997895148220428
pnorm(-1.96)         # 0.024997895148220428
pnorm(-1.96, false)  # 0.9750021048517796

for z = -25:0.1:25
    a = alnorm(z)
    p = pnorm(z)
    if abs(a - p) > 1e-9
        println("z = $z, a = $a, p = $p, diff = $(abs(a-p))")
    end
end

#=====
z = -1.2, a = 0.11506968113439031, p = 0.11506967022170828, diff = 1.0912682035790766e-8
z = -1.1, a = 0.13566606961261096, p = 0.13566606094638264, diff = 8.666228318299218e-9
z = -1.0, a = 0.15865526063064378, p = 0.15865525393145705, diff = 6.69918673312786e-9
z = -0.9, a = 0.18406013036234564, p = 0.1840601253467595, diff = 5.015586140855177e-9
z = -0.8, a = 0.21185540219053878, p = 0.21185539858339666, diff = 3.6071421127825687e-9
z = -0.7, a = 0.2419636546893566, p = 0.24196365222307298, diff = 2.466283621771481e-9
z = -0.6, a = 0.2742531193320773, p = 0.2742531177500736, diff = 1.582003694711176e-9
z = 0.6, a = 0.7257468806679227, p = 0.7257468822499265, diff = 1.5820037502223272e-9
z = 0.7, a = 0.7580363453106433, p = 0.758036347776927, diff = 2.466283621771481e-9
z = 0.8, a = 0.7881445978094612, p = 0.7881446014166034, diff = 3.6071421405381443e-9
z = 0.9, a = 0.8159398696376543, p = 0.8159398746532405, diff = 5.015586168610753e-9
z = 1.0, a = 0.8413447393693563, p = 0.8413447460685429, diff = 6.699186649861133e-9
z = 1.1, a = 0.864333930387389, p = 0.8643339390536173, diff = 8.666228290543643e-9
z = 1.2, a = 0.8849303188656097, p = 0.8849303297782918, diff = 1.0912682091301917e-8
=====#

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