#==========
Julia の修行をするときに,いろいろなプログラムを書き換えるのは有効な方法だ。
以下のプログラムを Julia に翻訳してみる。
母相関係数の差の検定に必要なサンプルサイズ・検出力
http://aoki2.si.gunma-u.ac.jp/R/power-cor-test.html
ファイル名: powercortest.jl 関数名: powercortest
翻訳するときに書いたメモ
cor0, cor1 を求めるときは,二通り求まる場合がある。
==========#
using Distributions, Roots
function powercortest(; n=NaN, cor0=NaN, cor1=NaN, siglevel=NaN, power=NaN,
alternative="two.sided")
sum(map(isnan, [n, cor0, cor1, siglevel, power])) == 1 ||
error("n, cor0, cor1, power, siglevel のどれか一つだけを NaN にする")
alternative in ["one.sided", "two.sided"] ||
error("'alternative' は 'one.sided' か 'two.sided'")
tside = alternative=="two.sided" ? 2 : 1
if isnan(power)
power = cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))
elseif isnan(n)
n = find_zero(n -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))-power, [4, 1e7], Bisection())
elseif isnan(cor0)
a = []
try
b = find_zero(cor0 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))-power, [-1, cor1], Bisection())
append!(a, b)
catch
println("No cor0 in [-1, cor1] can be found to achieve the desired power")
end
try
b = find_zero(cor0 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))-power, [cor1, 1], Bisection())
append!(a, b)
catch
println("No cor0 in [cor1, 1] can be found to achieve the desired power")
end
cor0 = length(a) == 1 ? a[1] : (a[1], a[2])
elseif isnan(cor1)
a = []
try
b = find_zero(cor1 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))-power, [-1, cor0], Bisection())
append!(a, b)
catch
println("No cor1 in [-1, cor0] can be found to achieve the desired power")
end
try
b = find_zero(cor1 -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))-power, [cor0, 1], Bisection())
append!(a, b)
catch
println("No cor0 in [cor0, 1] can be found to achieve the desired power")
end
cor1 = length(a) == 1 ? a[1] : (a[1], a[2])
elseif isnan(siglevel)
siglevel = find_zero(siglevel -> cdf(Normal(), sqrt(n-3)*abs(atanh(cor0)-atanh(cor1))-
cquantile(Normal(), siglevel/tside))-power, [1e-5, 0.99999], Bisection())
else
error("internal error")
end
println(" n = $n")
println(" cor0 = $cor0")
println(" cor1 = $cor1")
println(" sig. level = $siglevel")
println(" power = $power")
println("alternative = $alternative")
Dict(:n => n, :cor0 => cor0, :cor1 => cor1, :siglevel => siglevel,
:power => power, :alternative => alternative)
end
powercortest(siglevel=0.05, power=0.8, cor0=0.0, cor1=0.4) # find n
# n = 46.73160799446543
# cor0 = 0.0
# cor1 = 0.4
# sig. level = 0.05
# power = 0.8
# alternative = two.sided
powercortest(siglevel=0.05, power=0.8, cor0=0.3, cor1=0.4) # find n
# n = 605.5778584973955
# cor0 = 0.3
# cor1 = 0.4
# sig. level = 0.05
# power = 0.8
# alternative = two.sided
powercortest(n=605, siglevel=0.05, cor0=0.3, cor1=0.4) # find power
# n = 605
# cor0 = 0.3
# cor1 = 0.4
# sig. level = 0.05
# power = 0.7996236163502786
# alternative = two.sided
powercortest(n=605, power=0.8, cor0=0.3, cor1=0.4) # find siglevel
# n = 605
# cor0 = 0.3
# cor1 = 0.4
# sig. level = 0.050157266442876956
# power = 0.8
# alternative = two.sided
powercortest(n=605, power=0.8, siglevel=0.05, cor0=0.3) # find cor1
# n = 605
# cor0 = 0.3
# cor1 = (0.19288845143335967, 0.40004600000298063)
# sig. level = 0.05
# power = 0.8
# alternative = two.sided
powercortest(n=605, power=0.8, siglevel=0.05, cor1=0.4) # find cor0
# n = 605
# cor0 = (0.2999501647530865, 0.4913458910993233)
# cor1 = 0.4
# sig. level = 0.05
# power = 0.8
# alternative = two.sided