とにかく，なんでもかんでも，tidyverse のみを使って問題解決すべき（できないなら，あなたのスキルが足りないだけだ）というのは，新興宗教的（今ではコンプライアンス的に問題のある表現かもしれないが）
プログラミングには疎いが tidyverse なら何でも知ってるぞという人は別なんだけど。
を見ると，仰天する。
代替法の例が，以下に示されている。
よく，あんなたくさんの関数群を操れるなあ。
自分も，おなじような経験はある。自分の都合の良い関数を10数個作って，それらの組み合わせである程度のことをできるようにするシステムを作った。
数百個の関数を自由自在に組み合わせる方法を示して，「何でもできるよ！」というのは，本人以外にはほぼ理解不能ではないか?
逆に言えば，何でもできるように，新たな関数を次々に作る。
Julia Version 1.7.1 (2021-12-22) リーリースされた。
macOS ARM (M-series Processor) では，まだ RCall は使えない。
中澤先生が，
> RStudioはR-4.2.0 Prereleaseを呼び出そうとしてエラーが起こり，終了するしかなくなってしまうようだ
と書いていたので，そうだったかな?と思い，
colMeans(iris[1:4])
hist(iris[,1])
だけをやってみたけど，特に問題なく動いた。
RStudio 2021.09.0 Build 351
"Ghost Orchid" Release (077589bc, 2021-09-20) for macOS
Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
であった。かなり古いバージョンのようであるが，「No Update Available / You're using the newest version of RStudio」と言われる。
中澤先生は Windows でやったのかな?
多くの新機能，変更がある。
https://docs.julialang.org/en/v1/NEWS/
そのうちで，私が気になったものを以下に抜萃した。
√，∛ の前の乗算記号は不要になった
julia> x, y = 2, 3
(2, 3)
julia> x√y
3.4641016151377544
julia> x * sqrt(y)
3.4641016151377544
julia> x∛y
2.8844991406148166
julia> x * y^(1/3)
2.8844991406148166
julia> x * cbrt(y)
2.8844991406148166
julia> 3√3
5.196152422706632
julia> 3 * sqrt(3)
5.196152422706632
#####
配列の連結を指定するセミコロンの個数でどのように連結するかを指定できるようになった
julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
1 2
3 4
julia> B = [11 12; 13 14]
2×2 Matrix{Int64}:
11 12
13 14
julia> [A; B] # 1 個の場合は 1 次元方向へ
4×2 Matrix{Int64}:
1 2
3 4
11 12
13 14
julia> [A;; B] # 2 個の場合は 2 次元方向へ
2×4 Matrix{Int64}:
1 2 11 12
3 4 13 14
julia> [A;;; B] # 3 個の場合は 3 次元方向へ。以下同様に...
2×2×2 Array{Int64, 3}:
[:, :, 1] =
1 2
3 4
[:, :, 2] =
11 12
13 14
#####
文字列の途中でバックスラッシュ \ の後に改行し，次行以降に文字列を続ける場合，\ も改行コードも除去されて一連の文字列になる。次行以降の文字列にインデンテーションがある場合には画面上はインデンテーションされるが，文字列の中には空白もタブも含まれない。
julia> s = "abcdefg\
hijklmn"
"abcdefghijklmn"
julia> s = "abcdefg\
hijklmn"
"abcdefghijklmn"
#####
range() の引数指定法が改善された。
range(start, stop, length)
range(start, stop; length, step)
range(start; length, stop, step)
range(;start, length, stop, step)
#####
replace() で複数の置換ができるようになった。置換は左から右へ同時に行われる。
julia> replace("abcabc", "a" => "b", "b" => "c", r".+" => "a")
"bca"
julia> replace("abcabc", "a" => "b")
"bbcbbc"
julia> replace("abcabc", "a" => "b", "b" => "c")
"bccbcc"
#####
乱数発生器が Mersenne Twister から Xoshiro256++ になった。
このため，Ramdom.seed!() で乱数の種が設定されていても，Mersenne Twister のときと同じ乱数列は生成されないので注意が必要。
R, Python, Julia の速度比較 2021/12
実験環境
macOS Monterey バージョン 12.1
Mac mini (M1, 2020)
チップ Apple M1
メモリ 8GB
R M1 チップ対応
version 4.2.0
platform aarch64-apple-darwin21.1.0
arch aarch64
os darwin21.1.0
system aarch64, darwin21.1.0
Python M1 チップ対応
Python 3.10.0
Julia M1 チップ対応
Version 1.7.0 (2021-11-30)
使用するデータ
> n = 100000
> m = 1500
> x = data.frame(round(matrix(rnorm(n*m, 50, 10), n, m), 2))
> system.time({write.csv(x, file="test.csv", row.names=FALSE)})
user system elapsed
59.049 1.063 60.135
データフレームとして読み込み
R の場合
> system.time({df = read.csv("test.csv", colClasses="numeric")})
user system elapsed
12.105 0.906 13.627
ちなみに，data.table パッケージの fread だと更に速い。
> library(data.table)
> system.time({df = fread("test.csv")})
user system elapsed
1.723 0.258 3.058
しかし，fread で読み込まれるのは data.table なので，data.frame と若干の違いがあるので注意。
以下は，read.csv で読み込んだ data.frame を使う。
Python の場合
>>> from time import time
>>> import pandas as pd
>>> start = time(); df = pd.read_csv('test.csv'); print(time() - start)
6.904027938842773
Julia の場合
using CSV, DataFrames
@time df = CSV.read("test.csv", DataFrame); # 4.313486 seconds
@time df = CSV.read("test.csv", DataFrame, type=Float64); # 3.268237 seconds
単変量関数
合計
> system.time({s = colSums(df)})
user system elapsed
0.368 0.318 0.743
>>> start = time(); s = df.sum(); print(time() - start)
2.547268867492676
@time s = sum.(eachcol(df)); # 0.044548 seconds
算術平均
> system.time({m = colMeans(df)})
user system elapsed
0.353 0.084 0.436
>>> start = time(); m = df.mean(); print(time() - start)
0.2523488998413086
@time m = mean.(eachcol(df)); # 0.041935 seconds
不偏分散
> system.time({v = sapply(df, var)})
user system elapsed
0.487 0.002 0.491
>>> start = time(); v = df.var(ddof=1); print(time() - start)
0.6684608459472656
@time v = var.(eachcol(df)); # 0.080120 seconds
標準偏差
> system.time({sd = sapply(df, sd)})
user system elapsed
0.496 0.003 0.497
>>> start = time(); sd = df.std(ddof=1); print(time() - start)
0.7022347450256348
@time sd = std.(eachcol(df)); # 0.076600 seconds
中央値
median() はデータをソートする必要があるので，処理時間は長くなる。
> system.time({sapply(df, median)})
user system elapsed
2.086 0.364 2.691
>>> start = time(); s = df.median(); print(time() - start)
1.9600980281829834
@time med = median.(eachcol(df)); # 1.667501 seconds
二変量関数
ピアソンの積率相関係数
> system.time({cor(df[,1], df[,2])})
user system elapsed
0.005 0.001 0.005
>>> import numpy as np
>>> start = time(); r = np.corrcoef(df['X1'], df['X2'])[0,1]; print(time() - start)
0.0030248165130615234
>>> from scipy.stats import pearsonr
>>> start = time(); r = pearsonr(df['X1'], df['X2'])[0]; print(time() - start)
0.0011849403381347656
@time cor(df.X1, df.X2) # 0.000154 seconds
多変数関数
ピアソンの積率相関係数行列
> system.time({r = cor(df)})
user system elapsed
106.154 0.377 106.509
>>> start = time(); r = df.corr(); print(time() - start)
377.98001408576965
@time r = cor(Matrix(df)); # 3.985498 seconds
スピアマンの順位相関係数行列
> system.time({rs = cor(df[1:200], method="spearman")})
user system elapsed
4.698 0.152 4.858
>>> start = time(); rs = df.iloc[:, 0:200].corr(method='spearman'); print(time() - start)
23.894737243652344
@time rs = corspearman(Matrix(df[:, 1:200])); # 217.821656 seconds
ケンドールの順位相関係数行列
> system.time({rk = cor(df[1:5], method="kendall")})
user system elapsed
327.977 0.465 328.402
探索していて見つけたライブラリ
https://rdrr.io/cran/pcaPP/man/cor.fk.html
library(pcaPP)
> system.time({cor.fk(df[1:5])})
user system elapsed
0.106 0.002 0.109
> system.time({cor.fk(df[1:50])})
user system elapsed
10.330 0.262 10.560
>>> start = time(); rk = df.iloc[:, 0:5].corr(method='kendall'); print(time() - start)
0.2503321170806885
>>> start = time(); rk = df.iloc[:, 0:50].corr(method='kendall'); print(time() - start)
26.7141010761261
@time rk = corkendall(Matrix(df[:, 1:5])) # 0.134137 seconds
@time rk = corkendall(Matrix(df[:, 1:50])) # 11.786005 seconds;
分散共分散行列
> system.time({v = var(df)})
user system elapsed
106.169 0.688 106.949
>>> start = time(); v = df.cov(); print(time() - start)
3.9224328994750977
@time v = cov(Matrix(df)); # 6.487346 seconds
重回帰分析
> system.time({ans = lm(X1 ~ X2+X3+X4+X5+X6+X7+X8+X9+X10, data=df)
+ summary.ans = summary(ans)
+ coeff = ans$coeff
+ R2 = summary(ans)$r.squared})
user system elapsed
0.046 0.028 0.136
>>> from sklearn.linear_model import LinearRegression
>>> start = time();
>>> clf = LinearRegression()
>>> y = df['X1']
>>> x = df[['X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10']]
>>> ans = clf.fit(x, y)
>>> coef = ans.coef_
>>> intercept = ans.intercept_
>>> print(time() - start)
0.031842947006225586
using GLM
@time begin
result = lm(@formula(X1 ~ X2+X3+X4+X5+X6+X7+X8+X9+X10), df)
coefs = coef(result);
R2 = r2(result)
end # 0.025926 seconds
独立変数の個数が非常に多い場合
> system.time({
+ ans = lm(X1 ~ ., data=df)
+ summary.ans = summary(ans)
+ coeff = ans$coeff
+ R2 = summary(ans)$r.squared})
user system elapsed
168.687 2.262 172.028
>>> start = time();
>>> y = df['X1']
>>> x = df.drop('X1', axis=1)
>>> ans = clf.fit(x, y)
>>> coef = ans.coef_
>>> intercept = ans.intercept_
>>> R2 = ans.score(x, y)
>>> print(time() - start)
11.533896923065186
using GLM
@time fit(LinearModel, hcat(ones(size(df, 1)), Matrix(select(df, Not(:X1)))), df.X1); # 5.237764 seconds
JITコンパイラでPythonを高速化する「Pyjion」、バージョン1.0が登場
https://project.nikkeibp.co.jp/idg/atcl/19/00002/00281/
との記事があったが，それとは違うものかもしれないが，何もしなくても JIT コンパイラが動いているようなのだが。
以下の実行例で，np.mean() を 3 回実行しているが，初回は 0.0298 秒だが，2 回目以降は 0.0184 秒になっている。
Python 3.10.0 (v3.10.0:b494f5935c, Oct 4 2021, 14:59:19) [Clang 12.0.5 (clang-1205.0.22.11)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> np.random.seed(123)
>>> x = np.random.normal(50, 10, 100000000)
>>> from time import time
>>> start = time(); m = np.mean(x); print(m, time() - start)
50.000602597052335 0.029835939407348633
>>> start = time(); m = np.mean(x); print(m, time() - start)
50.000602597052335 0.01838517189025879
>>> start = time(); m = np.mean(x); print(m, time() - start)
50.000602597052335 0.018398046493530273
なお，x の長さを 10 倍にして， x = np.random.normal(50, 10, 1000000000) について同じことをやってみると，12.85 秒，14.40 秒，13.35 秒となった。これは，ガーベージコレクションやスワップが多発するために理論的には10倍になるが，遥かに経過時間が長くなっているのであろう。
ちなみに Julia も JIT コンパイラが働くが，
julia> using Random
julia> Random.seed!(123);
julia> @time x = randn(100000000);
0.241307 seconds (2 allocations: 762.939 MiB, 1.91% gc time)
julia> @time x = x .* 10 .+ 50;
0.099196 seconds (228.56 k allocations: 775.553 MiB, 2.98% gc time, 26.03% compilation time)
julia> @time mean(x)
0.052130 seconds (153.62 k allocations: 8.408 MiB, 46.63% compilation time)
50.001400590232656
julia> @time mean(x)
0.016738 seconds (1 allocation: 16 bytes)
50.001400590232656
julia> @time mean(x)
0.016801 seconds (1 allocation: 16 bytes)
50.001400590232656
であった。
長さを 10 倍したときにはそれぞれ，17.980272，13.516718，12.985337 秒であった。やはり，理論値よりは長い。
もう，Python なんて，いらね〜〜〜〜〜じゃない?
ここ数年?
満月のたびに，「今日は〇〇ムーン」とかネットで書かれてですね。
なんだそれ。
「アメリカインディアンはこの満月を〇〇と呼んでいました」
はあ，そうですか。私，日本人ですけど何か?
流星群とか，満月とか，一般人（失礼します）にアピールするのも良いのだけど。
例えば，流星群。期待するほどのものじゃないですよ。（まあ，「今まで，流れ星を見たことないんです❤」なんて人は別格なんですが。別格なら別格で，「なんだ，期待はずれじゃんプンプン」てこともあるので）
程々にしてほしいものです。
なお，夢を壊すかもしれないけど，散在流星というのもありましてですね，季節も時間も問わず，一時間に1つ流れるかどうかという流星もあるんです。これこそ，みられたら「幸せになれるかもれない」。流れ星が消えるまでに何回願い事を言えばいいのか?なんてこともない。それこそ，「まれなこと」なんだから。
夢が叶うかどうかと言うなら，流星群ではなく散在流星のほうがご利益ありそうですよね。
心を無にして，夜空を眺めましょう。
流れ星がみられなくても，あなたが見つめた煌めく星星は，あなたに流れ星以上のすばらしい何かを伝えてくれたのではないでしょうか
みんな大好き，Julia の速度自慢
https://blog.goo.ne.jp/r-de-r/e/35653aa771bcd15136f14dae7969e593
では触れなかったが，
みんな大好き，R と Python の速度比較
https://blog.goo.ne.jp/r-de-r/e/9d0ef78b17007c6de642983c368aaf20
で，R でのケンドールの順位相関係数の計算が遅すぎると書いた。
これを解決できるパッケージを見つけたので報告する。
実験環境
macOS Monterey バージョン 12.1
Mac mini (M1, 2020)
チップ Apple M1
メモリ 8GB
R M1 チップ対応
version 4.2.0
platform aarch64-apple-darwin21.1.0
arch aarch64
os darwin21.1.0
system aarch64, darwin21.1.0
Python M1 チップ対応
Python 3.10.0
Julia M1 チップ対応
Version 1.7.0 (2021-11-30)
テストに使用するデータ
R で作った，10万件1500変数のデータフレームを csv ファイルとして保存。
これを各プログラムでデータフレームとして入力して各種関数の実行速度を測定する。
n = 100000
m = 1500
x = data.frame(round(matrix(rnorm(n*m, 50, 10), n, m), 2))
write.csv(x, file="test.csv", row.names=FALSE)
df = read.csv("test.csv", colClasses="numeric")})
まずは，R
> system.time({rk = cor(df[1:5], method="kendall")})
user system elapsed
327.977 0.465 328.402
探索していて見つけたライブラリ
https://rdrr.io/cran/pcaPP/man/cor.fk.html
library(pcaPP)
> system.time({fk = cor.fk(df[1:5])})
user system elapsed
0.106 0.002 0.109
cor() より 3000倍速い。
> system.time({fk = cor.fk(df[1:50])})
user system elapsed
10.330 0.262 10.560
次いで， Python
>>> start = time(); rk = df.iloc[:, 0:5].corr(method='kendall'); print(time() - start)
0.2503321170806885
cor.fk より 2 倍以上遅い。
>>> start = time(); rk = df.iloc[:, 0:50].corr(method='kendall'); print(time() - start)
26.7141010761261
最後に，Julia
@time rk = corkendall(Matrix(df[:, 1:5])) # 0.134137 seconds
@time rk = corkendall(Matrix(df[:, 1:50])) # 11.786005 seconds;
cor.fk() より若干遅いが，Python に比べて 2倍以上速い。
curve(1/(1+exp(-x)), -4, 4)
using Plots
plot(x -> 1/(1+exp(-x)), label="")
かつては年金番号の統合ミスとか誤入力とか，今回はワクチン接種記録システムとか。
https://blog.goo.ne.jp/r-de-r/e/f8b88d3e3fa5f0a452bb0d39d2f6a854
https://blog.goo.ne.jp/r-de-r/e/ef134fa356cfdc59c9a39b06871614bc
https://blog.goo.ne.jp/r-de-r/e/ef134fa356cfdc59c9a39b06871614bc
https://blog.goo.ne.jp/r-de-r/e/f8b88d3e3fa5f0a452bb0d39d2f6a854
function pell(n)
x = isqrt(big(n))
(y, z, r, e1, e2, f1, f2) = (x, 1, 2x, 1, 0, 0, 1)
while true
y = r * z - y
z = (n - y * y) ÷ z
r = (x + y) ÷ z
(e1, e2) = (e2, e2 * r + e1)
(f1, f2) = (f2, f2 * r + f1)
(a, b) = (f2, e2)
(b, a) = (a, a * x + b)
a * a - n * b * b == 1 && return (a, b)
end
end
n = 61
a = pell(n) # (1766319049, 226153980)
a[1]^2 - n * a[2]^2 # 1
n = 109
a = pell(n) # (158070671986249, 15140424455100)
a[1]^2 - n * a[2]^2 # 1
n = 181
a = pell(n) # (2469645423824185801, 183567298683461940)
a[1]^2 - n * a[2]^2 # 1
n = 277
a = pell(n) # (159150073798980475849, 9562401173878027020)
a[1]^2 - n * a[2]^2 # 1
n x n 魔方陣で，n が奇数の場合の生成法
https://ja.wikipedia.org/wiki/魔方陣#バシェー方式
using Printf
function f(n)
a = zeros(Int, n, n)
m = n ÷ 2
for i = 1:n
(sx, sy) = (i - 1 - m, n - i - m)
for j = 1:n
(x, y) = (sx + j, sy + j)
x <= 0 && (x += n)
x > n && (x -= n)
y <= 0 && (y += n)
y > n && (y -= n)
a[x, y] = (i - 1) * n + j
end
end
for i = 1:n # 以下は結果の表示
for j = 1:n
@printf("%5i", a[i, j])
end
println()
end
a
end
a = f(15)
106 219 92 205 78 191 64 177 50 163 36 149 22 135 8
9 107 220 93 206 79 192 65 178 51 164 37 150 23 121
122 10 108 221 94 207 80 193 66 179 52 165 38 136 24
25 123 11 109 222 95 208 81 194 67 180 53 151 39 137
138 26 124 12 110 223 96 209 82 195 68 166 54 152 40
41 139 27 125 13 111 224 97 210 83 181 69 167 55 153
154 42 140 28 126 14 112 225 98 196 84 182 70 168 56
57 155 43 141 29 127 15 113 211 99 197 85 183 71 169
170 58 156 44 142 30 128 1 114 212 100 198 86 184 72
73 171 59 157 45 143 16 129 2 115 213 101 199 87 185
186 74 172 60 158 31 144 17 130 3 116 214 102 200 88
89 187 75 173 46 159 32 145 18 131 4 117 215 103 201
202 90 188 61 174 47 160 33 146 19 132 5 118 216 104
105 203 76 189 62 175 48 161 34 147 20 133 6 119 217
218 91 204 77 190 63 176 49 162 35 148 21 134 7 120
sum(a, dims=1)
1×15 Matrix{Int64}:
1695 1695 1695 1695 1695 1695 … 1695 1695 1695 1695 1695 1695
sum(a, dims=2)'
1×15 adjoint(::Matrix{Int64}) with eltype Int64:
1695 1695 1695 1695 1695 1695 … 1695 1695 1695 1695 1695 1695
using LinearAlgebra
sum(tr(a)) # 1695
sum(tr(rotr90(a))) # 1695
我慢のならないことを，我慢ならん!といい続けたいものだ。たとえ，それが多くの人に受け入れられなくても。
基礎を大切に。
