Julia の変数スコープに注意

R の fisher.test() において，hybrid=true を選択した場合に使用される rcont2() を Julia に移植した。

Julia は，for ブロック，while ブロックが 新たな変数スコープを持つので，注意が必要な局面がある。

普通は，エラーメッセージで気づくのだが，REPL 環境でやっているとおかしなことになってしまい，ドツボに嵌まることもある。

この関数でも，3 箇所（3 変数）で注意が必要。初期化は必要ない（どのような数値を代入しても無関係）だが，そこでその変数を定義しないと「そんな変数ないよ」エラーになる。試行錯誤・テストのために PEPL で適当な値を設定してしまっていれば，このエラーが出なくなってします。そこで，ドツボにはまる（for や while の外に出たとたんに変数の値が代わってしまうという現象が起きる）。

もうひとつは，Julia らしい，面白いプログラミング技法

・ 論理式 && 論理式が真の場合に実行される文

・ 論理式 || 論理式が偽の場合に実行される文

if a < 5
     break
end

なんてのが，一行で書ける（わかりにくくなってしまうかもしれないが，慣れると便利）

using SpecialFunctions

function rcont2(nrowt, ncolt)
  #=============================================================================
  RCONT2 constructs a random two-way contingency table with given sums.
  Licensing:
    This code is distributed under the GNU LGPL license.
  Modified:
    28 Janualy 2021
  Author:
    Original FORTRAN77 version by WM Patefield.
  Reference:
    WM Patefield,
    Algorithm AS 159:
    An Efficient Method of Generating RXC Tables with
    Given Row and Column Totals,
    Applied Statistics,
    Volume 30, Number 1, 1981, pages 91-97.
  Parameters:
    Input, Array{Int64,1}nrowt(nrow), Array{Int64}ncolt(ncol),
      the row and column sum vectors.
    Output, Array{Int64,2}matrix(nrow, ncol), the matrix.
  =============================================================================#

  MAXTOT = 1000000
  nrow = length(nrowt)
  ncol = length(ncolt)
  ntotal = sum(ncolt)
  nrow > 1 || error("The length of the row sum vector must be greater than 1.")
  ncol > 1 || error("The length of the column sum vector must be greater than 1.")
  all(nrowt .> 0) || error("An entry in the row sum vector must be greater than 0.")
  any(ncolt .> 0) || error("An entry in the column sum vector must be greater than 0..")
  ntotal == sum(nrowt) || error("The row and column sum vectors must have the same sum.")
  MAXTOT >= ntotal || error("The sum of the column sum vector entries is too large.")

  # Calculate log-factorials.
  fact = logfactorial.(0:ntotal)
  # Construct a random matrix.
  matrix = zeros(Int64, nrow, ncol)
  jwork = zeros(Int64, ncol)
  jwork[1:ncol-1] = ncolt[1:ncol-1]
  jc = ntotal
  ib = 99999 # 初期化は不要だが，定義は必要
  for l = 1 : nrow - 1
    nrowtl = nrowt[l]
    ia = nrowtl
    ic = jc
    jc -= nrowtl
    for m = 1 : ncol - 1
      id = jwork[m]
      ie = ic
      ic -= id
      ib = ie - ia
      ii = ib - id
      # Test for zero entries in matrix.
      if ie == 0
        ia = 0
        matrix[l, m:ncol] = 0
        break
      end
      # Generate a pseudo-random number.
      r = rand(1)[1]
      # Compute the conditional expected value of MATRIX(L,M).
      done1 = false
      nlm = 99999 # 初期化は不要だが，定義は必要
      while true
        nlm = floor(Int64, ia * id / ie + 0.5)
        iap = ia + 1
        idp = id + 1
        igp = idp - nlm
        ihp = iap - nlm
        nlmp = nlm + 1
        iip = ii + nlmp
        x = exp(fact[iap] + fact[ib+1] + fact[ic+1] + fact[idp] -
          fact[ie+1] - fact[nlmp] - fact[igp] - fact[ihp] - fact[iip])
        r <= x && break
        sumprb = x
        y = x
        nll = nlm
        lsp = false
        lsm = false
        # Increment entry in row L, column M.
        done2 =false # 初期化は不要だが，定義は必要
        while ! lsp
          j = (id - nlm) * (ia - nlm)
          if j == 0
            lsp = true
          else
            nlm += 1
            x *= j / (nlm * (ii + nlm))
            sumprb = sumprb + x
            if r <= sumprb
              done1 = true
              break
            end
          end
          done2 = false
          while !lsm
            # Decrement the entry in row L, column M.
            j = nll * (ii + nll)
            if j == 0
              lsm = true
              break
            end
            nll -= 1
            y = y * j / ((id - nll) * (ia - nll))
            sumprb = sumprb + y
            if r <= sumprb
              nlm = nll
              done2 = true
              break
            end
            !lsp && break
          end
          done2 && break
        end
        (done1 || done2) && break
        r = sumprb * rand(1)[1]
      end
      matrix[l, m] = nlm
      ia -= nlm
      jwork[m] -= nlm
    end
    matrix[l, ncol] = ia
  end
  # Compute the last row.
  matrix[nrow, 1:ncol-1] = jwork[1:ncol-1]
  matrix[nrow, ncol] = ib - matrix[nrow, ncol-1]
  return matrix
end

# 使い方

rcont2([1,2,3], [2,4])