統計解析R たぶんpart3くらい
ググって出てきた最初のコード試して見たけどJIS丸めと同じ結果だね。 結局きちんと動作するコードは見つからなかったよ 少し驚いたな 四捨五入は簡単なようでいて落とし穴のある難問なのかもね 0.285 には以下の通り2進数に起因する誤差があるから、 内部的には 0.284999..... 扱いなんだろうね > 29 - (0.285*100 + 0.5) [1] 3.552714e-15 確認には直接フォーマット系の関数を使った方が簡単だと思う。 ``` {r} sprintf ('%.64f', 0.15); ``` 多くの四捨五入の実装は`10`を掛ける演算を使っているが、注意が必要になる。 C99の標準関数を使えば、`round (x, digit = 0)`に相当する四捨五入は得られる。 ``` {Rcpp} #include <Rcpp.h> #include <cmath> // [[Rcpp::plugins(cpp17)]] // [[Rcpp::export]] double std_round (double x) { return std::round (x); } ``` `std_round`を使って`round (x, digit = 1)`に相当することをやろうとして、 次のようなコードを書くと、オワコンがやって来る。 ``` {r} sprintf ('%.64f', std_round (0.15 * 10) / 10); ``` 数値計算と表示との整合性がとれていないにも関わらず、通常のプリント文を 見ている限り気が付きにくい厄介なバグになる。敗因は次の計算にある。 ``` {r} sprintf ('%.64f', 0.15 * 10); ``` 四捨五入のスマートな実装方法は思いつかない >>364 >>366 内部レジスタの値は正しく丸められてても 表示でこけてるだけって可能性もあるからな 疑うのは大事だけど 自分が勘違いしてないことを確認することも大事 >>370 デフォルト環境で0.29が0.28と表示されるわけないわ 以下のページでRの丸めについて英語で説明してるね ttps://cran.r-project.org/web/packages/round/vignettes/Rounding.html 最終的に多倍長演算ライブラリ gmp の使用を推奨してる library(gmp) myround <- function(x, digits){ return(as.numeric((x*10^digits+0.5)/(10^digits))) } val <- as.bigq(285, 1000) # 0.285 を 285/1000 と分数指定 myround(val, digits=2) そういえばPython2と3でなんかそんな問題があったような・・・ 2ではできたのに3ではできなくなったみたいな。 >>372 myround に floor (負の扱いによっては trunc) 付け忘れてた 補完して解釈してください sprintfでの確認やgmpライブラリの方法は勉強になるなあ gmpを使った四捨五入は完璧! ただbigqに変換するところの手作業が実用上のネックで、これを無くすのが難しそう 残念ながら多倍長数への変換の自動化は無理そうだ 手作業で一つひとつ変換するとなると直接手作業で四捨五入したほうがずっと早いしミスも少ない Rで実用的な四捨五入は不可能なのだろうか 数値が "±nn.nnn" 形式の文字列だけなら以下の関数でどうだろ ("2.85e-1" とかだと駄目だけど) library(gmp) text2bigq <- function(txt){ denominator <- 10^nchar(txt) # 10^小数点以下の桁数より大きい整数 numerator <- as.integer(round(as.numeric(txt)*denominator)) return(as.bigq(numerator, denominator)) } myround <- function(x, digits){ return(as.numeric(floor(x*10^digits+0.5)/(10^digits))) } val <- text2bigq("0.285") myround(val, digits=2) text2bigqなんてよく考えつくなあ これ引数のクオーテーション無くてもいけるね それなら2.85e-1でも通る ただ数値によって桁あふれのエラーが出るのが難点だね 丸め桁に対して数値の桁数が十分多いときは分岐して普通にroundすればいけるのかな? ついにR4.1.0来たね 目玉はネイティブパイプの実装 iris |> head() 積と[作用](https://ncatlab.org/nlab/show/action )の関係 ``` {r} self = with (new.env (), { id = function (a) a; mult = function (bc) function (ab) function (a) bc (ab (a)); flip = function (abc) function (b) function (a) abc (a) (b); r_curry_fw = function (abc) function (a) function (b) abc (a, b); r_curry_bw = function (abc) function (a, b) abc (a) (b); `%>%` = r_curry_bw (flip (id)); `%*%` = r_curry_bw (flip (mult)); f = function (x) function (y) paste0 ('(', x, ', ', y, ')', sep = ''); testthat::test_that ('', { testthat::expect_equal ( 1 %>% f (2) %>% f (3) # 作用 , 1 %>% (f (2) %*% f (3)) # 積 ); }); as.list (rlang::current_env ()); }); ``` ユーザー定義の二項演算をサポートしないプログラミング言語では、 パイプは言語仕様に組み込まないと実現できないので、パイプを言語仕様に 取り入れるか否かという議論がぼちぼちある。その一方で、関数の合成については 眼中無しになっている。積を積極的にサポートするプログラミング言語が 少ないのは何故だろう? Rに限った話ではないので、何らかの人間の特性が関わっているのかもしれない。 上のコードの空白や改行も含めた全ての文字中で、約20パーセントを"function" という単語が占める。今度導入されたラムダ表記だと、バックスラッシュの 占める割合は約3パーセントにまで減る。個人的には、これが嬉しい。 パイプのオマケでの仕様変更な気がするが、塞翁が馬。 ネイティブパイプは実際に使ってみると思ってたよりいまいちだった うーん cars |> head() |> {\(x) lm(dist ~ speed, x)}() cars |> head() |> x => lm(dist ~ speed, x) やはりこっちが楽 cars %>% head %>% lm(dist ~ speed, .) ネイティブパイプって任意の引数に値を渡せず第一引数限定なんだ。 それは昔ながらの引数が統一されていない関数だとちょっと使いにくいな。 これでもできたけどやはりめんどくさい かなり見やすくはなった cars |> head() |> lm(formula = dist ~ speed) Rとは関係しないかもしれない: 1. [JavaScript and the next decade of data programming (2020) | Hacker News](https://news.ycombinator.com/item?id=27373388 ) インタラクティブな可視化 1. [Learn R Through Examples (2020) | Hacker News](https://news.ycombinator.com/item?id=27404103 ) 賛否両論 1. [Yann LeCun Deep Learning Course 2021 | Hacker News](https://news.ycombinator.com/item?id=27387154 ) ありがたや 1. [246B, Notes 3: Elliptic functions and modular forms | What's new](https://terrytao.wordpress.com/2021/02/02/246b-notes-3-elliptic-functions-and-modular-forms/ ) 四捨五入も"周期的"な関数`round (x + n) = round (x) + n` ``` {python round_py} import numpy; import mpmath; def round_theta (x, h = 1e-6): q = mpmath.exp (- 0.5 / h); z = lambda x: x / 2j / h; x = [mpmath.re ( mpmath.jtheta (3, z (x), q, 1) / 2j / mpmath.jtheta (3, z (x), q, 0) ) for x in x]; return numpy.array (x, dtype = numpy.float64); ``` インタラクティブな可視化やアニメーションのプロットをやってみたいけど印刷してハンコが前提なせいで機会がない インタラクティブな可視化はどうしてもGUIとの連携が必要になるが、 アニメは連携が要らない。 ``` {r, animation.hook = "gifski"} self = with (new.env (), { plot_text = \(...) \(text) { c01 = c (0, 1); plot (c01, c01, type = 'n', axes = F, ann = F); text (mean (c01), mean (c01), text, ...); }; purrr::map (strsplit ("あいうえお", "") [[1]], plot_text (cex = 3)); plot_text (cex = 5) ("\U1F44C"); as.list (rlang::current_env ()); }); ``` `gifski`ならアニメ専用のコードを必要としない。このことは一長一短だと思うが、 手っ取り早い。 functionふたつ繋げてmapを受けやすくしてるのか なるほどそういう方法もあるんだ わかっていないところなので最初に結論: * 混ぜるな危険 * 触らぬ神に祟りなし [NSE](http://adv-r.had.co.nz/Computing-on-the-language.html )を使った関数と [高階関数](https://en.wikipedia.org/wiki/Higher-order_function )を混ぜるな。 殆どの人にとって、何行でプログラムが書けるかは問題でなく、何分で結果が 得られるかが問題だと思う。見つけにくいバグを生みやすいテクは避けた方が 無難だと思う。 複数の表に対して同じ操作をするオモチャを考える。 ``` {r sub, dependson = ""} sub = with (new.env (), { r_swap = \(abc) \(b, a) abc (a, b); try_it = \(expr) tryCatch (expr, error = gettext); as.list (rlang::current_env ()); }); ``` つづく ``` {r map, dependson = "sub"} sub = with (sub, { `%map%` = purrr::map; f = \(abc) \(...) \(a) abc (a, ...); g = \(abc) \(b) \(a) abc (a, b); datum = list (head (cars), tail (cars)); datum %map% f (`[`) (2 : 1) |> try_it () |> print (); datum %map% g (`[`) (2 : 1) |> try_it () |> print (); datum %map% f (r_swap (lm)) (dist ~ speed) |> try_it () |> print (); datum %map% g (r_swap (lm)) (dist ~ speed) |> try_it () |> print (); datum %map% f (subset) (dist == max (dist)) |> try_it () |> print (); datum %map% g (subset) (dist == max (dist)) |> try_it () |> print (); datum %map% f (dplyr::group_by) (dist) |> try_it () |> print (); datum %map% g (dplyr::group_by) (dist) |> try_it () |> print (); datum %map% f (dplyr::mutate) (b = dist + 1) |> try_it () |> print (); datum %map% g (dplyr::mutate) (b = dist + 1) |> try_it () |> print (); as.list (rlang::current_env ()); }); ``` * `f (group_by) (dist)`と`g (group_by) (dist)`の違いは? * `lm`と他のNSEとの違いは? 個人的には、これらの疑問を明朗会計できない。 おしまい グループ化はg()だけエラーになるのか 意識すらしてなかった部分だけど面白い検討結果だね つまるところtidyverseの処理はtidyverse的に書くのが短くて早くてわかりやすいってことかな 4.1になってWindowsでもtidyxlでxlsxを読めるようになっていることに気付いた tidyxlはあまり使わないものの散らかったxlsx読むときには便利 前説が長いので最初に結論: >353は真似しないで欲しい。 [IHaskell](https://github.com/gibiansky/IHaskell )でコードを書く。 ``` haskell import GHC.Exts (groupWith) 分割 = groupWith length 仕事 (a : as) = (length a, length as + 1) 集計 = foldl push mempty where push out a = out ++ pure a 個々 = map 物件 = words "Lorem ipsum dolor sit amet, consectetur adipiscing elit" 比較 前処理 後処理 = lhs == rhs where lhs = 後処理 $ 集計 . 個々 仕事 . 分割 <$> 前処理 物件 rhs = 後処理 $ 集計 <$> 個々 仕事 <$> 分割 <$> 前処理 物件 -- 比較 Identity runIdentity 比較 (\a -> [a, a ++ a]) (flip (>>=) id) 比較 (\a j -> [a !! k | k <- j]) (flip id [1, 3, 5]) ``` [関手](https://en.wikipedia.org/wiki/Functor )と呼ばれるデザインパターンを 使っている。このパターンをRに翻訳する。 ``` {r self, dependson = ""} self = with (new.env (), { id = \(a) a; hom_snd = \(bc) \(ab) \(a) bc (ab (a)); flip = \(abc) \(b) \(a) abc (a) (b); const = \(a) \(b) a; r_curry_fw = \(abc) \(a) \(b) abc (a, b); r_curry_bw = \(abc) \(a, b) abc (a) (b); `%.%` = r_curry_bw (hom_snd); as.list (rlang::current_env ()); }); ``` つづく 例題を次の記事から拝借する。[fmap]{#fmap} * [Split-Apply-Combine and Map-Reduce in R](https://burtmonroe.github.io/SoDA501/Materials/SplitApplyCombine_R/ ) ``` {r fmap, dependson = "self"} self = with (self, { 分割 = \(x) x |> split (~ cyl); 仕事 = \(x) coef (lm (mpg ~ wt, x)); 集計 = \(x) purrr::map_dbl (x, \(x) x ["wt"]); 個々 = \(f) \(x) purrr::map (x, f); 物件 = mtcars; 比較 = \(管, 前処理, 後処理) with (list (`%管%` = 管), { lhs = 前処理 (物件) %管% (集計 %.% 個々 (仕事) %.% 分割) |> 後処理 (); rhs = 前処理 (物件) %管% 分割 %管% 個々 (仕事) %管% 集計 |> 後処理 (); testthat::test_that ("", testthat::expect_equal (lhs, rhs)); # lookme a = 前処理 (物件); b = 管 (a, 分割); c = 管 (b, 個々 (仕事)); d = 管 (c, 集計); rhs = 後処理 (d); testthat::test_that ("", testthat::expect_equal (lhs, rhs)); }); 比較 (r_curry_bw (flip (id)), \(x) rbind (head (x), tail (x)), id); 比較 (purrr::map, \(x) list (head (x), tail (x)), id); 比較 (r_curry_bw (flip (hom_snd)), \(x) \(y) rbind (head (x), y), \(x) x (tail (物件))); as.list (rlang::current_env ()); }); ``` つづく 1つめの"比較"は表を面倒くさい方法で`lm`している。2つめの"比較"は 表の代わりに表のリストを渡している。この場合の"%管%"が前回の投稿の `%pipe%`に相当する。3つめの"比較"は表の代わりに表への写像を渡している。 このように関手パターンを使うと、引数が様々な"形状"に変化する。 本題に入る。以前、次のようなコードを書いた。[boxplot]{#boxplot} ``` {r, dependson = ""} none = with (new.env (), { . = ggplot2::ggplot (mtcars); . = . + ggplot2::aes (x = factor (cyl), y = wt); . = . + ggplot2::geom_boxplot (); print (.); }); ``` `.`が定数である限り問題はないが、途中一箇所でも`. `が関数になると、 無限ループする。[fmap](#fmap)で`lookme`以下の変数`a..d`を全て`a`にすると、 三番目の"比較"で無限ループする。厄介なことに、全て`a`にしても、 一番目と二番目の"比較"は普通に動作する。関数内の変数は検索経路が 異なることに起因する。 単純化すると、次のコードは動くが、 ``` {r, dependson = ""} tryCatch ({ a = 1; a = 2 + a; a; }, error = identity) |> print (); ``` つづく 次のコードは無限ループする。 ``` {r, dependson = ""} tryCatch ({ a = 1; a = \(.) . + a; a (2); }, error = identity) |> print (); ``` PythonとJSも似たような挙動をするので、現在主流のインタープリターでは [boxplot](#boxplot)のような横着コードはご法度かもしれない。 おしまい 特典がいろいろありそう。 * [Modern Text Features in R](https://www.tidyverse.org/blog/2021/02/modern-text-features/ ) 話変わって>384の続き まずコードを並べる。 ``` {r dict, dependson = "round_py"} dict = with (new.env (), { who_max_h = \(h) \(x) { x = x - max (x); x = exp (x / h); x / sum (x); }; who_min_h = \(h) \(x) who_max_h (h) (- x); ge_h = \(h) \(x) 0.5 * (1 + tanh (0.5 * x / h)); le_h = \(h) \(x) ge_h (h) (- x); as.list (rlang::current_env ()); }); ``` つづく ``` {r round_r, dependson = "dict"} dict = with (dict, { round_theta = \(h) \(x) { reticulate::py $ round_theta (reticulate::np_array (x), h); }; round_1nn = \(h, n) { self = (- n) : n; \(x) purrr::map_dbl (x, \(x) { x = self - x; p = who_min_h (h) (0.5 * x * x); sum (self * p); }); }; round_count = \(h, n) { self = 1 : n - 0.5; \(x) purrr::map_dbl (x, \(x) { sum (ge_h (h) (x - self) - ge_h (h) (- x - self)); }); }; as.list (rlang::current_env ()); }); ``` つづく ``` {r round_anime, dependson = "round_r", animation.hook = 'gifski', fig.width = 4 * 3} dict = with (dict, { x_max = 3; x = x_max * seq (- 1, 1, len = 1e+2); x_huge = x_max * seq (- 1, 1, len = 1e+3); h = 10 ^ c (- 8, - 6, - 4, (- 3) : 3); old_par = par (mfrow = c (1, 3)); on.exit (par (old_par)); purrr::map (h, \(h) { if (h < 1) { x = x_huge; } title = sprintf ('h = %.1e', h); draw = \(y, ylab, ...) { plot (x, y, type = 'l', ylim = c (- x_max, x_max), xlab = 'x', ylab = ylab, main = title, ...); }; draw (round_theta (h) (x), 'theta'); draw (round_1nn (h, x_max) (x), 'deep'); draw (round_count (h, x_max) (x), 'shallow'); }); as.list (rlang::current_env ()); }); ``` つづく 書けるかな? 次の関数から始める。 ~~~ {.r} nearest_neighbor_integer = \(n) { self = (- n) : n; \(x) purrr::map_dbl (x, \(x) { x = self - x; j = which.min (0.5 * x * x); self [j]; }); }; ~~~ 実数`x`に最も近い`(-n):n`中の整数を選んでいる。この関数の`which.min` のところを`(- n):n`上の確率分布に変更したものを`round_1nn (h, n) (x)` としている。`n`無限大の極限をとると、 [ヤコビのテータ関数](https://en.wikipedia.org/wiki/Theta_function ) で書けて`round_theta (x, h)`になる。ヤコビの三重積を有限和で近似して `round_count (h, n) (x)`を得る。`n`を"容量"、`h`を"温度"と書く。 `ge_h (h) (x)`は低温極限でデジタル的な関数`ifelse (x >= 0, 1, 0)`になり、 `round_count`の低温極限は、容量が十分大きければ、ゼロと`x`の間にある 半整数`n + 1 / 2`の数を数える関数になる。 [1-nn](https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm ) で解いていた問題`round_1nn`が、ヤコビの三重積を通して、数え上げの問題 `round_count`に転化した形になっている。 アニメを見ると、低温では全部同じ感じだが、高温で`round_theta`とその他に 違いが出てくる。対称性`round_theta (x + k, h) == round_theta (x, h) + k` にガードされて、`round_theta`だけは、高温極限でもゼロへの定数写像に ならない。他の2つは、容量を有限で近似しているために、この周期性が 成り立たない。 つづく 四捨五入や浮動小数点の話を蒸し返す感じになるけど… 1.275-1と0.275を第二位でround()してみる 前者は0.27で後者は0.28になる x1 <- c(1.275-1, 0.275) round(x1, 2) [1] 0.27 0.28 sprintf()で眺めてみるとなんとなく理由がわかる気がする sprintf("%.24f", x1) [1] "0.274999999999999911182158" "0.275000000000000022204460" 0.265でも試してみると両者とも0.26になる x2 <- c(1.265-1, 0.265) round(x2, 2) [1] 0.26 0.26 sprintf()で眺めてみると… sprintf("%.24f", x2) [1] "0.264999999999999902300374" "0.265000000000000013322676" 偶数丸めで切り下げになるケースだけ誤差の範囲が広いのかな?? 書けるかな?書けたら話をぶった切ってごめん。who_max_hとge_hは個別に 覚えるより、次の関数から芋づる式に覚えた方が安上がりかもしれない。 ニョロニョロニョロ ドット ハスケル max_h h a b = h * log (exp (a / h) + exp (b / h)) min_h h a b = - max_h h (- a) (- b) ニョロニョロニョロ max_hの低温極限はRのpmaxになる。関数 reduce_max (h, x) := reduce (.init = - Inf, .x = x, .f = max_h)は log_sum_expと呼ばれ、低温極限がRのmaxになる。reduce (0, x, +)にsumを 使うのと同じで、"モノイドはreduceしとけ"パターンの1つになっている。 reduce_max (h, x)をxについて微分すると、soft_maxになり、低温極限がRの which.maxをワンホットで表したものになる。コードではsoft_maxをwho_max_h と書いている。関数max_zero (h, x) := max_h (h, x, 0)の低温極限はreluの 有限温度版にあたるsoft_plusと呼ばれる。max_zero (h, x)をxについて微分 すると、シグモイド関数になる。コードではシグモイド関数をge_hと書いている。 ge_h (h, x)をxについて微分すると、Rのdlogisになる。どの関数も低温極限で デジタル的な関数になり、アルゴリズム的な描像を与える。 温度hをプランク定数だと思うと、マスロフの脱量子化という言い方も理解できる。 シグモイド関数はフェルミ分布関数とも呼ばれるが、対となるボース分布関数は zero_maxの逆写像を微分して得られる。ボース分布関数を使っても、 ヤコビの三重積経由で、また別の低温極限が四捨五入になる関数が作れる。 reduce_maxの連続版はlog_integral_expとなるが、reduce_minの連続版は ラプラス近似と呼ばれる。ラプラス近似の補正第一項を正規分布にとることが 多いので、max_zeroが正規分布に対応すると見ることもできる。また、 max_zero (h, reduce_max (h, x))をxについてルジャンドル変換すると、 soft_maxの負のエントロピーになる。 おしまい どうせ過疎だから話をぶった切るなんて気にせず好きに書き込んでいいと思うよ 最近知ったpurrrとpatchworkの組み合わせは役に立たないけど面白いと思った library(patchwork) library(tidyverse) iris %>% nest_by(Species) %>% mutate(plot = list(ggplot(data, aes(Sepal.Length, Petal.Length)) + geom_point() + ggtitle(Species))) %>% pull(plot) %>% reduce(`+`) https://i.imgur.com/Qet93Lj.png ルジャンドル変換を一般化した凸共役についてまとまった表がある。 * [凸共役](https://en.wikipedia.org/wiki/Convex_conjugate ) この記事中の表の下二行と次の記事を照らし合わせてみる。 1. [Softmax Cross-Entropy and Negative Sampling](https://arxiv.org/abs/2106.07250 ) 2. [Noise Contrastive Estimation](https://arxiv.org/abs/1202.3727 ) | 凸 | 凸共役 | Bregman情報量 | |:--|:--|:--| | ln (1 + e^x) | x^* ln x^* + (1 - x^*) ln (1 - x^*) | KL情報量 | | - ln (1 - e^x) | ln p (x^*) + x^* ln (1 - p (x^*)) | NCE | p (x) := 1 / (1 + x) 前に書いた四捨五入と対比すると、KLの方がフェルミ分布関数、NCEの方が ボーズ分布関数に対応する。四捨五入の場合はフェルミ分布関数とボーズ分布関数 の裏にテータ関数がいた。KLとNCEの裏にも何かいるかな? 論文をちゃんと読めば答えもしくはヒントがあるかもしれない。 JIS丸めの話がどうしてこうなった?JIS丸め恐るべし。 JIS丸めや四捨五入については皮肉なことにExcelのほうがずっと信用できたりする 良くも悪くもうまいこと誤差処理してくれてる 402の数値をExcelで四捨五入するとすべて期待どおりの結果が得られる ggrepelってのいいね 凡例と見比べなくてすむ図が作れる これならカテゴリ多くて色が飽和する場合でも対処できそう 極端な例 https://i.imgur.com/3aKWSES.png readr2.0.0から複数ファイルの読み込みが可能になっていた map_dfr()使わなくていいし読み込みも圧倒的に速くて良い readr 2.0.0のwrite系関数は、カラム名に日本語があると正しく処理されないバグがないか? write_csv試したけど列名に日本語あっても問題なかった でもデータも日本語にしたらたしかに化けた ちなWindows すべての言語は、Linux 用 例外は、Windows, Android, iOS 用の言語だけ >>414 追確認ありがとう。やっぱり化けるか。 カラム名の出力抑止したら日本語含むデータは無事に出力されたので、当面、これで逃げとく。 >>416 write_csvのヘルプ見たらutf8で保存されるとあったからこれが理由ぽい データと列名の両方をutf8にしたら化けずに出力できた 具体的には以下を噛ませた #データ変換 %>% mutate(across(where(is.character), enc2utf8)) %>% #列名変換 %>% rename_with(enc2utf8) %>% RStudioチートシート一式が更新されてる dplyrなんかは構文が少し違ってきてたからありがたい >>413 自己レス readr 2.0.1にアップデートしたら>>417 の変換を行わないでもヘッダー、データに日本語を含む場合でも意図した通りに出力されるようになった。 >>419 これ見て自分も2.0.1にした 確かに問題なくなった 報告ありがとう 間違って機械学習スレに書き込んでしまったので、こちらに再書き込み。 Rmdファイル名に日本語が入っているとRStudioでread系の出力が表示されなくなるという謎のトラップに嵌った。 read.csv関数、readr::read_csv関数、readxl::read_excel関数などを端折ってread系と書いてしまった。 動くんだけどエラーも何も表示されないという現象だった。 適当に日本語のrmdからread_csv試してみたけど動作もメッセージも特に問題ないみたい 自分が条件を勘違いしてるだけかもだけど readrの2.1.0が来てた デフォルトで遅延読込するのが評判悪くて戻したらしい Rstudioで100個のtxtファイル (中身はTsvで3万行2列、1列目は全ファイル共通) を3万行、101列のマトリクスファイルにしようとしています。 left_joinで1列目をキーにして結合するためにデータを100個読み込む際に lf <- list.files(full.names = T) data <- lapply(lf,read.delim) を用いたのですが、 mat <- full_join(data)で結合させようとすると 'full_join' をクラス "list" のオブジェクトに適用できるようなメソッドがありません となってしまいました。 別の関数でこのような結合が可能でしょうか。 >>426 dplyr::bind_cols() があるとわかり(rowsしか見つけられていませんでした)、 こちらでdata2<-dplyr::bind_cols(data[1],data[2]) でエラーが出ないことまで分かったのですが、 今度は 1列目が全残りしてしまいました。。。 >>426 詳細はわかりませんがとりあえず読み込みの部分を以下にすれば結合されたdata.frame(tibble)ができると思います data <- readr::read_tsv(lf) あ、ごめんなさい 勘違いしてたのでいったん428は無視してください bind_colの代わりにdplyr::left_join()で結合してください >>426 確認だけど各ファイルの一列目は言わるインデックスで、どのファイルも全く同じ値で全く同じ並びなの? 例えば、1,2,3,...30,000のような。 >>430 コメントいただきありがとうございます。 data2<-dplyr::left_join(data[1],data[2]) UseMethod("left_join") でエラー: 'left_join' をクラス "list" のオブジェクトに適用できるようなメソッドがありません となってしまいました。 >>431 そのご認識であっています >>432 1列目がすべて同じ値でデータを結合させるキーとして使う必要がなければ、こんな感じ。 最後に1列目をバインドかジョインさせれば、概ね目的が達成できるのではないかと。 require(tidyverse) # サンプルデータの作成 for (i in 1:10) { data.frame( a = seq(from = 1, to = 100), b = rnorm(100) ) %>% readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) } # サンプルデータの結合 list.files("./sample_data", full.names = TRUE) %>% purrr::map_dfc(.f = function(x){readr::read_csv(x) %>% dplyr::select(2)}) 列名でデータが識別出来るようにしたければ、もうひと工夫必要だけど。 ファイル名を列の識別に使うとこんな感じ # サンプルデータの読み込みと結合 list.files("./sample_data", full.names = TRUE) %>% purrr::map_df(.f = function(x){readr::read_csv(x) %>% dplyr::mutate(file = basename(x))}) %>% tidyr::pivot_wider(names_from = file, values_from = b) >>433 おお!ありがとうございます。 解釈しながら進めてみます。 >>434 の方は、1列目をキーにleft_joinするのと同じ結果になるので、こちらの方が処理としては汎用性が高いかも。 # サンプルデータの作成 for (i in 0:4) { data.frame( a = seq(from = 1, to = 100), b = rnorm(100) ) %>% readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) } # サンプルデータの作成(aの値の範囲を少し変えてある) for (i in 5:9) { data.frame( a = seq(from = 51, to = 150), b = rnorm(100) ) %>% readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) } # サンプルデータの読み込みと結合(aが一致しない場合は欠損値NAとなる) list.files("./sample_data", full.names = TRUE) %>% purrr::map_df(.f = function(x){readr::read_csv(x) %>% dplyr::mutate(file = basename(x))}) %>% tidyr::pivot_wider(names_from = file, values_from = b) >>427 一列目がすべて同じとのことなのであれば、427の続きで data2[ , -seq(3,199,2)] とすれば必要な列だけ抽出できると思います あとはマトリクスに変換してください join系でうまく行ってなかったのはリストの扱いがわからなかったからなのですね リストにアクセスする場合、 data[1]でなくdata[[1]]としてデータフレームを取り出します >>437 すみません。本題ではないと思うのですが、 > for (i in 1:10) { + data.frame( + a = seq(from = 1, to = 100), + b = rnorm(100) + ) %>% + readr::write_excel_csv(file = paste0("./sample_data/sample_", i, ".csv")) のところで エラー: Cannot open file for writing:* './sample_data/sample_1.csv' が出てしまいました・・・ >>438 先にディレクトリ(フォルダ)作らないとだめよ。 なお、そのコードはサンプル用のデータファイルを作るだけだから、既にデータファイルがあるなら作らなくても構いません。 >>438 繰り返しになりますが、427にあるようにbind_cols(data)でdata2まで作成できたのであれば、あとは437の式で重複するindex列を消せば抽出が完了します マトリクス形式にするにはas.matrix()を使います エラーが出たという433のコードについては439さんがもう回答してくれていますが念のため補足説明します そのコードは作業ディレクトリにsample_dataという名前のフォルダを作ってから試す必要があります Windowsならエクスプローラからフォルダを作成すれば良いでしょう Rでコードを実行するとそのフォルダの中にテスト用のcsvが作成されます(フォルダを覗いてみてください) そのテスト用csvを使った汎用的な方法が色々と紹介されていますので、どういう挙動をするか試してみてください あとは自分の実データに応用すればOKです cranにR 4.2.0がきてた Windowsもutf8に マルチバイトの苦労が減るかな 利点欠点は知らないけどとりあえず試しに移ってみては Vscodeに移行完了。 他の言語と環境同じなのは楽。簡単に複数のプロセスを走らせることができて便利。 環境変数とかの設定は大変。 RStudioもPositだかに名称変更して多言語対応を前面に出すらしいから色々迷うところだね tidyverseにlubridateが追加されたらしい これは便利かも ALTREPを勉強したノート: https://www.klgrth.io/paste/n3aj8 有効期限の選択肢が最長2日で、2日経つと消えちゃう。 一応HTMLということでアップしたが、HTMLとして表示できない。 ダウンロードしてブラウザーで開くと、多分、表示できると思う。 中身はRというより、殆どC++とJSで、偶に日本語が入るという感じになっている。 facet_wrapで強制的に任意の行×列にする方法がわからない ggh4xのfacet_wrap2はx軸が別の値に変わる致命的なバグがあって駄目だった cowplotとかで強引にやるしかないかな こういう事? [r - What's the difference between facet_wrap() and facet_grid() in ggplot2? - Stack Overflow](https://stackoverflow.com/questions/20457905/whats-the-difference-between-facet-wrap-and-facet-grid-in-ggplot2 ) ```{r} g <- ggplot2::ggplot (ggplot2::mpg, ggplot2::aes (displ, hwy)) g + ggplot2::facet_grid (cyl ~ class) g + ggplot2::facet_wrap (cyl ~ class) ``` 前回のALTREPから離れて配列型のイテレーターに進んでみた。 https://www.klgrth.io/paste/wcdat 今回のファイル形式はRmdにしてみた。 Rcppが使える環境なら動作するんじゃないかと思う 例えばmpgデータならcylでもclassでも常に3*3のサイズにfacetしたい(不足分はスペース) 結局patchworkのplot_layout()が良かった facet_wrap2()よりかなりめんどくさいが仕方ない RのヒープからC++のコンテナーで使うメモリを確保するようにした: https://www.klgrth.io/paste/ampar 今回はサードパーティのライブラリを使っているので、そのままでは コンパイルできないと思う。 * [GitHub - martinus/unordered_dense: A fast & densely stored hashmap and hashset based on robin-hood backward shift deletion](https://github.com/martinus/unordered_dense ) ヘッダーファイル一枚のライブラリなので、Rmdの中にヘッダーファイルを コピペすればコンパイルできるようになると思う。 前回のカウンターの続きで最終回 [zwzg9](https://www.klgrth.io/paste/zwzg9 ) 今回は`ankerl`に加えて`boost`も使っている。`boost`はヘッダー一枚 というわけにはいかないので、面倒かもしれない。使っているのは、 `hash_value`と`hash_combine`という2つの関数だけなので、適当な関数で 差し替えられると思う。`NA`と`NaN`の関係を調べている節は、 intelエンディアンのみの対応で、armエンディアンでは動作しないと思う。 実を言うと、`NaN`がいっぱいあることを知らなかった。1954は何の年だろう? ALTREPを勉強した時にかなりギットハブを徘徊した。そこが勉強のピークだっと 思う。今回の実装では、殆どRの勉強はせずに、知ってる関数だけを使っている。 1. `TYPEOF / ALTREP` 1. `Rf_xlength` 1. `DATAPTR_OR_NULL / DATAPTR / XXX_ELT` 1. `Rf_allocVector` 違うアプリケーション、例えば、評価系だと違う関数セットになって、 また勉強が必要になるかもね。 ファセットの問題は{facetious}のfacet_wrap_strict()でもできた こっちは問題なさそう [facet_wrap_strict](https://github.com/coolbutuseless/facetious )と `facet_wrap`を並べてみる。 ``` {r na_1954, dependson = ""} . = dplyr::mutate (mtcars, cyl = factor (cyl, lev = c (1 : 8, 1000L))) . = ggplot2::ggplot (., ggplot2::aes (mpg, wt)) . = . + ggplot2::geom_point () . + ggplot2::facet_wrap (~ cyl, nrow = 3, ncol = 3) . + facetious::facet_wrap_strict (~ cyl, nrow = 3, ncol = 3) . + ggplot2::facet_wrap (~ cyl, nrow = 3, ncol = 3, drop = F) ``` 神エクセル繋がりということで [A Bayesian probability worksheet](https://terrytao.wordpress.com/2022/10/07/a-bayesian-probability-worksheet/ ) その例で行くなら… cylはいじらずそのままでfacet_wrap_strictするのが求める図かな RのS5クラス使って継承する場合、コンストラクタ中で継承するクラスのコンストラクタを呼ぶことはできないのかな? つまりクラスAとクラスBを定義して、クラスBがクラスAを継承している場合、クラスBのコンストラクタ中でクラスAのコンストラクタを呼ぶ方法 ちょっとググったらS4クラスは出来そうなんだが、S5クラスのは解決策が見当たらんかった 別途自前でクラスAのコンストラクタ本体を定義して、クラスBのコンストラクタからそれを呼ぶしかないか?! dplyrを1.1.0にしてみた .byはまあまあ便利だけどtally()には使えなかった ナンノコッチャと思って調べたら [A Heisenbug lurking in async Python | Hacker News](https://news.ycombinator.com/item?id=34754276 ) 専門用語だったでござる。 [Heisenbug - Wikipedia](https://en.wikipedia.org/wiki/Heisenbug ) read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる