Julia に翻訳--199 因子分析の適合度検定

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

因子分析の適合度検定
http://aoki2.si.gunma-u.ac.jp/R/fa.fit.test.html

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

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

計算の途中で結果が Matrix{Any} になることがある。今回の場合は

sc = [i == j ? 1 / sqrt(uniquenesses[i]) : 0 for i = 1:p, j = 1:p]

では，sc が Matrix{Real} になるのが原因。それを回避するには

sc = [i == j ? 1.0 / sqrt(uniquenesses[i]) : 0.0 for i = 1:p, j = 1:p]

とすればよいようだが，なんだかなあ！

==========#

using DataFrames, Statistics, StatsBase, Rmath, LinearAlgebra

function fafittest(data, factors, uniquenesses)
    n, p = size(data)
    r = cor(data)
    sc = [i == j ? 1.0 / sqrt(uniquenesses[i]) : 0.0 for i = 1:p, j = 1:p]
    r = sc * r * sc
    values = eigvals(r) # Julia ではデフォルトで昇順
    e = values[1:p - factors]
    s = -sum(log.(e) .- e) + factors - p
    chisq = (n - 1 - (2 * p + 5) / 6 - (2 * factors) / 3) * s
    df = ( (p - factors) ^ 2 - p - factors) / 2
    pvalue = pchisq(chisq, df, false)
    println("chisq. = $chisq,  df = $df,  p value = $pvalue")
    Dict(:chisq => chisq, :df => df, :pvalue => pvalue)
end

data = [
  935  955  926  585 1010  925 1028  807  769   767
  817  905  901  632 1004  950  957  844  781   738
  768  825  859  662  893  900  981  759  868   732
  869  915  856  448  867  874  884  802  804   857
  787  878  880  592  871  874  884  781  782   807
  738  848  850  569  814  950  957  700  870   764
  763  862  839  658  887  900 1005  604  709   753
  795  890  841  670  853  874  859  701  680   772
  903  877  919  460  818  849  884  700  718   716
  761  765  881  485  846  900  981  728  781   714
  747  792  800  564  796  849  932  682  746   767
  771  802  840  609  824  874  859  668  704   710
  785  878  805  527  911  680  884  728  709   747
  657  773  820  612  810  849  909  698  746   771
  696  785  791  578  774  725  932  765  706   795
  724  785  870  509  746  849  807  763  724   760
  712  829  838  516  875  725  807  754  762   585
  756  863  815  474  873  725  957  624  655   620
  622  759  786  619  820  769  807  673  698   695
  668  753  751  551  834  849  807  601  655   642];

uniquenesses = [0.005, 0.186849474704386, 0.235609633298087,
  0.594548313805094, 0.005, 0.005, 0.644059344661751,
  0.005, 0.466720965448361, 0.573892982555005]
fafittest(data, 4, uniquenesses)
# chisq. = 15.163572526523414,  df = 11.0,  p value = 0.1751303147464357