Julia でスピアマンの順位相関係数行列

ご存じのこととは思いますが，念の為

Julia で 2 変数間のスピアマンの順位相関係数を求めるには，

using StatsBase, Statistics
corspearman([1,2,3,4,5], [3,2,1,2,3])

です。

スピアマンの順位相関係数行列を求めるとき，

function spearmancor0(x::Array{Float64,2})
    nc = size(x, 2)
    r = ones(nc, nc)
    for i = 2:nc
        for j = 1:i
            r[i, j] = r[j, i] = corspearman(x[:, i], x[:, j])
        end
    end
    return r
end

としてはいけません。

corspearman は データ行列も受け付けます。
以下のように定義した関数がやるように，データ行列を列ごとに平均順位を付け，そのあとで cor 関数を呼びます。

function spearmancor(x::Array{Float64,2})
    rank = similar(x, Float64)
    for i = 1:size(x, 2)
        rank[:, i] = tiedrank(x[:, i])
    end
    return cor(rank)
end

のようにしましょう。

時間計測すると，以下のようになります。

x = randn(10000, 100);
@time b = spearmancor0(x); # 8.367428 seconds (70.69 k allocations: 2.260 GiB, ...
@time a = spearmancor(x);  # 0.097368 seconds (714 allocations: 38.250 MiB, ...
@time c = corspearman(x); # 0.094456 seconds (1.06 k allocations: 38.337 MiB
a == b # true
a == c # true

90 倍ほど速いです。当然ですが，得られる相関係数行列は同じです。

セイウチ演算子(:=)って今更...（その2）

いやあ，うっかりしてたわ。

R でセイウチ演算子での代入箇所に複文が書けるって書いていたのに，Julia のこと書かなかったわ。

Julia では，複文は begin 〜 end で囲むので次の様にして，複文を書くことができる。

どんなに複雑な式を何個書いても問題ない。

が，もっとも大事な注意があることを書いてなかったわ。

a = [1,2,3,4,5]
if begin n = length(a); w=20*pi end > 10
    println("List is too long ($n elements, expected <= 10)")
else
    println("variable n is $n, w = $w")
end
println(n, ", ", w)

List is too long (5 elements, expected <= 10), w = 62.83185307179586
5, 62.83185307179586

おや，5 elements, expected <= 10 になってますよ。

あわてちゃいけない。複文の評価は，最後の文の評価結果なの！！つまり， w > 10 をやったことになってるの。

位置を考えるか，最後に n を置くか。以下のようにする。

if begin n = length(a); w=20*pi; n end > 10

＝＝＝＝＝＝＝＝＝＝

R も同じ。n > 10 をやりたいなら，複文の最後に n を置く。

a = c(1,2,3,4,5)
if ({n = length(a); w = 20*pi; n} > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10), w = %f\n", n, w))
} else {
    cat(sprintf("variable n is %d, w = %f\n", n, w))
}
cat(n, w)

variable n is 5, w = 62.831853
> cat(n, w)
5 62.83185

 

セイウチ演算子(:=)って今更...

Python 3.8 になって，使えるようになって感涙にむせんでいる人もいるようなんだけど...

a = [1,2,3,4,5]
if (n := len(a)) > 10:
    print(f"List is too long ({n} elements, expected <= 10)")
else:
    print(f"variable n is {n}")
print(n)

実行結果

variable n is 5
5

しかし，これと同じ機能は多くの言語では古くから使えていた。
というか，使えていたので，使い方を間違って問題を起こすこともあった。

R では，:= を使うときにも = を使えば良い。ただし，'(代入文)' のように '( )' でくくることを忘れてはいけない。

a = c(1,2,3,4,5)
if ((n = length(a)) > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10)\n", n))
} else {
    cat(sprintf("variable n is %d\n", n))
}
print(n)

実行結果

variable n is 5
> print(n)
[1] 5

R では if 文の条件式を ( ) でくくらなければならないので，Python と同じようには見えないかも知れないが，外側の( ) を取り除けば Python のセイウチ演算子と同じなのだ。

Julia だともっとわかりやすいだろう（条件式を ( ) でくくる必要がない）。
Python のセイウチ演算子の位置には，代入演算子 = を使っているだけだ。
セイウチ演算子なんてなくても，ちっとも困らない。

a = [1,2,3,4,5]
if (n = length(a)) > 10
    println("List is too long ($n elements, expected <= 10)")
else
    println("variable n is $n")
end
println(n)

実行結果

variable n is 5
5

古いというほど古くはない C でも以下のようになる。
代入式を(n = sizeof(a)/sizeof(int)) のように ( ) でくくるだけ。

#include
#include
int main(void) {
      int n;
      int a[5] = {1,2,3,4,5};
      if ((n = sizeof(a)/sizeof(int)) > 10) {
            printf("List is too long (%d elements, expected <= 10)\n", n);
      }
      else {
            printf("variable n is %d\n", n);
      }
      printf("n = %d\n", n);
      return EXIT_SUCCESS;
}

実行結果

variable n is 5
n = 5

ところで，R はもっと進んでいるんだ。

a = c(1,2,3,4,5)
if ({n = length(a)} > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10)\n", n))
} else {
    cat(sprintf("variable n is %d\n", n))
}
print(n)

のように，代入文を { } でくくっても良い。
それは，この位置に複文を書けるということだ。

以下のようなこともできる。すきなだけ複数の複雑な式を書ける。

a = c(1,2,3,4,5)
if ({n = length(a); m = 2*pi} > 10) {
    cat(sprintf("List is too long (%d elements, expected <= 10)\n", n))
} else {
    cat(sprintf("variable n is %d\n", n))
    cat(sprintf("2*pi = %f\n", m))
}
cat(n, m)

実行結果

variable n is 5
2*pi = 6.283185
> cat(n, m)
5 6.283185

しかし，同じように { } で複文が書ける C では，この位置に複文を書くとエラーになる。

もっと古い AWK でも同じように1つの文は書けるが，やはり複文は書けない。

BEGIN {
      a = 2
      b = 6
      if ((n = a * b) > 0) {
            print "true", n
      }
      else {
            print "false", n
      }
      print "outer", n
}

実行結果

true 12
outer 12

要するに，新しい演算子を作るかどうかということではなく，文法がどうか(記述されたものをどのように解釈するか）ということになるだろう。

数学の問題を Julia で解く

某所に Python による解答例があったので，Julia でやるならどうなるか書いてみた。

Q. 1.　2018n ≡ 2 mod(1000) を満たす正の整数nの最小値を求めよ。

for n = 1:99999
  if mod(2018n, 1000) == 2
    println(n) # 389
    break
  end
end

Q. 2.　tan(2 Arctan (1/3)+Arctan(1/12))を求めよ。

println(tan(2atan(1//3)+atan(1//12))) # 0.8888888888888888

a = tan(2atan(1//3)+atan(1//12))
println(Int(10a - a) // 9) # 8//9

Q. 3.　xyz空間内の4点(1,-4,1),(2,2,2),(2,-6,-3),(3,-2,-1)を頂点に持つ四面体の体積を求めよ。

sarasu(x1, y1, z1, x2, y2, z2, x3, y3, z3) =
    abs((y1*z2*x3+z1*x2*y3+x1*y2*z3-z1*y2*x3-x1*z2*y3-y1*x2*z3) / 6)

a1 = [1, -4, 1];
b1 = [2, 2, 2];
c1 = [2, -6, -3];
d1 = [3, -2, -1];
x1,　y1,　z1,　x2,　y2,　z2,　x3,　y3,　z3 = vcat(b1 .- a1, c1 .- a1, d1 .- a1);
println(abs((y1*z2*x3+z1*x2*y3+x1*y2*z3-z1*y2*x3-x1*z2*y3-y1*x2*z3) / 6)) # 3.0
println(sarasu(x1, y1, z1, x2, y2, z2, x3, y3, z3)) # 3.0

別解

A = hcat(b1 .- a1, c1 .- a1, d1 .- a1);
using LinearAlgebra
println(abs(det(A)) / 6) # 3.0

Q. 4.1　3次正方行列 A=[4 1 1; 1 2 1; 1 1 2] について固有値を求めよ。

A = [4 1 1; 1 2 1; 1 1 2];
using LinearAlgebra
println(round.(Int, eigvals(A))) # [1, 2, 5]

Q. 4.2　E=[1 0 0; 0 1 0; 0 0 1] のとき，A^3-8A^2+18A-12E を求めよ

E = Array{Int}(I, 3, 3);
println(A^3 .- 8A^2 .+ 18A .- 12E) # [2 1 1; 1 0 1; 1 1 0]

Q. 5　2つの機械A, Bが、1年間に故障する回数は，平均がそれぞれ2.5回、4.5回のポアソン分布に従う。A,Bの故障は独立である。
Q. 5.1　機械Aが1年間に4回以上故障する確率を求めよ。

using Distributions
println(round(ccdf(Poisson(2.5), 3), digits=3)) # 0.242

Q. 5.2　機械A,Bが合わせて8回以上故障する確率を求めよ。

A = Poisson(2.5)
B = Poisson(4.5)
answer = 1.0
for i = 0:7
    for j = 0:i
        answer -= pdf(A, i-j) * pdf(B, j)
    end
end
println(round(answer, digits=3)) # 0.401