統計解析R たぶんpart3くらい
[ブログ](https://colinfay.me/js-const-r/ )の [続き](https://twitter.com/_ColinFay/status/1176425719653642242 ) 必要があれば[R Online](https://srv.colinfay.me/r-online )で試せる ストーカーの有無によっても変わる ``` {r} with (list (), { a <- list (x = 1); lockBinding ('a', environment ()); tryCatch (a $ x <- 999, error = print); a; }); with (list (), { a <- list (x = 1); lockBinding ('a', environment ()); stalker <- a; tryCatch (a $ x <- 999, error = print); a; }); ``` https://twitter.com/5chan_nel (5ch newer account) Rでリストをつくってデータフレームを組むと X <- c(A, B, C) Y <- c(1, 2, 3) df <- data.Frame(X, Y) df A 1 B 2 C 3 って出てくると思ったら ABC123 って出てきます。これはどうしたらいいのでしょうか? これで、ちゃんとデータフレームになるけど? X <- c("A", "B", "C") Y <- c(1, 2, 3) df <- data.frame(X, Y) df X Y 1 A 1 2 B 2 3 C 3 [mathjax-node-page](https://github.com/pkra/mathjax-node-page ) を使って`knitr`が作ったMathJaxの数式をsvgに変換する 工具ディレクトリ`$(tool)`に`mathjax-node-page`をインストールする ~~~ {.bash} cd $(tool) npm install mathjax-node-page ~~~ ディレクトリ`$(tool)/node_modules/`ができる ファイル`$(tool)/page.js`を作る ~~~ {.js} const fs = require ('fs'); const path = require ('path'); const mathjax = require ('mathjax-node-page'); const argv = process.argv; const usage = () => { return ['[usage]', path.basename (argv [0]), path.basename (argv [1]), '<input-html-file>'].join (' '); }; if (argv.length !== 3) { console.log (usage ()); process.exit (1); }; const encoding = {encoding: 'utf-8'}; const inn_file = argv [2]; const inn = fs.readFileSync (inn_file, encoding); mathjax.mjpage (inn, {format: ['TeX']}, {svg: true}, (out) => { process.stdout.write (out); }); ~~~ よっこらしょ ~~~ {.bash} nodejs $(tool)/page.js R/a.html > html/a.html ~~~ オフラインでも読めるようになる もっと賢い方法があるかも Electron = Chromium(ブラウザ) + Node.js(サーバー・ローカルPC にアクセス) + V8(JavaScript エンジン) Electronが想像どおりのものだとすると アプリケーションを作る際には良い選択肢の一つだと思う ただ上でやろうとしていることに対しては真逆な方向に進んでしまうと思う HTMLをオフラインで読めるようにするにはRmdのYAMLを設定すれば良い [HTML document](https://bookdown.org/yihui/rmarkdown/html-document.html ) ~~~ {.yaml} output: html_document: mathjax: local self_contained: false ~~~ この方法で作ったHTMLをディレクトリごと配布すれば問題ない しかし完結した一枚のHTMLにまとめられればそれに越したことはない このちっさい野望を叶えるために Pandocは`self-contained`というオプションを用意している この行く手を拒むのがMathJaxになっている この話題は誰も興味がないと思うのでもう少し続ける Pandocのフィルターを使う手もある [mathjax-pandoc-filter](https://github.com/lierdakil/mathjax-pandoc-filter ) ~~~ {.yaml} output: html_document: pandoc_args: [--filter, $(tool)/node_modules/.bin/mathjax-pandoc-filter] ~~~ フィルターが確実に動作する限りはベストだと思う 多分`mathjax-node-page`より実行速度が早い 多分ブラウザ内でMathJaxを動作させるより実行速度が早い ただしフィルターにバグがあってPandocが解釈できないjsonを返すと HTMLが出力されない ちなみに最近のPandocではipynbも出力できる 次のように書いておけばipynbのコードブロックに変換してくれる ~~~ code (lambda who: print ('hello ' + who)) ('Jupyter'); ~~~ 逆にipynbをマークダウンに変換もしてくれる 使い方はいろいろだと思う おしまい 整然データって実験条件とかはどのように格納するべきですか? 1. 実験条件ファイル expDescr101.csv 氏名 性別 年齢 室温 モニタ解像度 問題数 問題パターン ・・・ TaroSat M 32 32 1920 29 A ・・・ 2. 実験データファイル 時刻 反応A 反応B 1011 1 0 1020 0 1 1100 0 0 以下続く 1の実験条件ファイルの条件が50個ぐらいあり、横にすごく長くなってしまうんですが、 1のファイル形式としてはこれで適切なのでしょうか? それだけじゃなんとも言えんが、ぱっと見て条件データと実験データってどうやって紐付けんの?とは思った。 >>323 1. 〇〇Descr.csv 2. 〇〇Data.csv というファイル名にしてファイル名で紐付けしようと思ってました。 Rだと1,2は同一ファイルに記載するものなのでしょうか? エクセルとかなら横の数気にするけど、 Rになるとどうでもいい感じw >>324 後からどうとでもなるけど、同一の方が楽だろうね。 検査者間信頼性としてICCを求めたいのですが、 ICC(1,1) (1,k)〜(3,k)まで6この下位モデルすべて求められる一般的なパッケージってありませんか? 調べたところ改変Rコマンダーってのがあるんですけど、 これって論文とかに使って問題ないでしょうか? ubuntu20.04のxfce環境だとRstudioのプルダウンメニューがマウスで開けない・・・・・ 仕方ないからキーボードでやってるけど同様の症状の人いる? Windowsだとファイル名のパスが通らないライブラリが多くて困る collapseとかいうライブラリーすごく速くていいな 以前書いたプルルのパイプのバグ(?)はフィックスしている。 [Tidyverse](https://www.tidyverse.org/blog/2020/11/magrittr-2-0-is-here/ ) いつフィックしたのかはわからないが、magrittrを2.0にアップデートしたら 動くようになっていた。次の記事を見て、新しく導入されるっぽい関数宣言の 短縮形を試していた際に気がついた。 [R adds native pipe and lambda syntax](https://news.ycombinator.com/item?id=25316608 ) Rockerさまさま。というわけでJSネタ。 [A Modern JavaScript Tutorial](https://news.ycombinator.com/item?id=25333350 ) JSのプロミス(Rのプロミスとは別)を勉強した時に、素人目線の注意書きが ところどころにあって、とても助かった記憶がある。皆さんを心地よい睡眠へ 誘う渾身のプレゼン資料を作る時に役立つかもしれない。 [[r-hub/r-minimal: Minimal Docker images for R](https://github.com/r-hub/r-minimal ) `rocker`のイメージとサイズを比べてみる。 ~~~ {.bash} REPOSITORY TAG SIZE rhub/r-minimal devel 35.3MB rocker/r-ver devel 813MB rocker/tidyverse devel 2.31GB ~~~ `rhub/r-minimal`にRのライブラリをインストールしようとすると、 `gcc`などのシステムのライブラリが必要になる。それを調べて インストールする手間を考えると、初めから全込みのイメージを ダウンロードした方が手っ取り早いかもしれない。それでも、 `r-minimal`から始めたくなる場面がしばしばあると思う。特に、 初めてDockerを使う場合は、まずは、`r-minimal`で試してみる価値が あると思う。金は時なり。 マイクロソフトがMicrosoft365Rというパッケージを公開してる 365(Office365)をRから操作できるらしい rvestのバージョンがついに1.0.0になったみたいよ 追加されたhtml_text2()はどんななんだろ タイトル勝ちと思ったけど... [The Hintons in your Neural Network](https://www.reddit.com/r/MachineLearning/comments/m2qo5p ) 物理では粒子をトン付で呼ぶ習慣がある。エレクトロン、フォトン、等。 で、[ヒントン](https://en.wikipedia.org/wiki/Geoffrey_Hinton ) redditではこの渾身のオヤジギャグが通じていないっぽい。 物理現象だと因果がはっきりしてるから○○エフェクトってなりやすいんだろうね 心理の場合でもプラシーボ効果やカリギュラ効果みたいに因果が明確なものは○○効果だね バタフライ効果はバタフライエフェクトと言われがちなのはなんなのだろう RStudioで普通に日本語入力できるようになってるんだね 以前は変換確定するまで数文字しか表示されなかったのに 下表のようにある職種米にあqるskillに対する得点を計測し平均値を表にしました。 この得点からスキルをABCの3つにクラスタリングしています。 医師 看護 介護 クラスタ skill1 0.8 0.4 0.2 A skill2 0.3 0.5 0.1 B skill3 0.2 0.3 03 C skill4 0.1 0.3 03 C ・・・ skill15 0.1 0.2 0.8 B で、クラスタAに対する医師・看護・介護の得点をボックスプロットしたいのですが、 なにかいい方法ありませんでしょうか? 下表のようにある職種米にあqるskillに対する得点を計測し平均値を表にしました。 この得点からスキルをABCの3つにクラスタリングしています。 医師 看護 介護 クラスタ skill1 0.8 0.4 0.2 A skill2 0.3 0.5 0.1 B skill3 0.2 0.3 03 C skill4 0.1 0.3 03 C ・・・ skill15 0.1 0.2 0.8 B で、クラスタAに対する医師・看護・介護の得点をボックスプロットしたいのですが、 なにかいい方法ありませんでしょうか? 宿題かなんか? 医師・看護・介護が水準になるtidy data形式にしてから描けばいいんじゃ? >>348 tidyverseは使えますか? データフレーム名をdfとすれば大枠は以下でいけますよ library(tidyverse) df %>% filter(クラスタ == "A") %>% pivot_longer(-クラスタ) %>% ggplot(aes(name, value)) + geom_boxplot() 上の回答が答えになっていると思う。以下はおまけ。 ``` {r nasdaq} df = " CEO, COO, CTO, cluster 0.8, 0.4, 0.2, A 0.3, 0.5, 0.1, B 0.2, 0.3, 03, C 0.1, 0.3, 03, C 0.1, 0.2, 0.8, B "; df = read.csv (text = df, colClasses = c ("numeric", "numeric", "numeric", "character"), strip.white = T); print (df); ``` tidy無し版 ``` {r plain, dependson = 'nasdaq'} for (clazz in unique (df $ cluster)) { go = subset (df, cluster == clazz, select = - cluster); print (go); boxplot (go, xlab = 'role', ylab = 'skill', main = sprintf ('cluster = %s', clazz)); } ``` 続く パイプ無し版 ``` {r tidy, dependson = 'nasdaq'} for_plot = tidyr::pivot_longer (df, - cluster, names_to = "role", values_to = "skill"); print (for_plot); for (clazz in unique (for_plot $ cluster)) { go = dplyr::filter (for_plot, cluster == clazz); print (go); . = ggplot2::ggplot (go); . = . + ggplot2::aes (x = role, y = skill); . = . + ggplot2::geom_boxplot (); . = . + ggplot2::labs ( title = sprintf ('cluster = %s', clazz) ); print (.); } ``` プロットを並べる。 ``` {r arrange_plot, dependson = 'tidy', fig.width = 4 * 2, fig.height = 4} . = ggplot2::ggplot (for_plot); . = . + ggplot2::aes (x = role, y = skill); . = . + ggplot2::geom_boxplot (); . = . + ggplot2::facet_wrap (~ cluster); print (.); ``` 箱を並べる。 ``` {r arrange_box, dependson = 'tidy', fig.width = 4 * 2, fig.height = 4} . = ggplot2::ggplot (for_plot); . = . + ggplot2::aes (x = role, y = skill, fill = cluster); . = . + ggplot2::geom_boxplot (); print (.); ``` おしまい 色々とご親切にありがとうございます 自分のつたいないコードよりスッキリ問題解決しました m(_ _)m 白状すると、検索してコピペして動いたコードを貼っただけ。 `boxplot`自体使ったことがなかった。 箱ひげ図は義務教育で教えるようになってからじわじわ市民権を得つつある感じ 10年前はほとんど誰も知らなくて見せてもなんだコレだったけど RでJIS丸めじゃない四捨五入ってどうやればいいの? round(0.285,2)が0.28じゃなく0.29になるような 検索して出てきたコードをいくつか試したけど2進数のせいか どれもround2(0.285,2)が0.28になってしまう このぐらいの桁数の丸めは需要ないのかな ググって出てきた最初のコード試して見たけど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系関数は、カラム名に日本語があると正しく処理されないバグがないか? read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる