統計解析R たぶんpart3くらい
データをこねくり回してグラフ書くにはexcelより融通が利いてしかも早い ような気がする >>175 そりゃそうでしょうね Excelはこねくり回さないで使うためのツールだから グラフはggplot2が強力だよね facetとかはExcelだと大変 デフォルトのテーマが少し古くさいけど 機械学習とかビッグデータとかPythonと何かの組み合わせに圧倒されててRとかどマイナーになってしもうた 適材適所でしょうね。 元々設計自体がデカいメモリーがあったあらと言う前提で 設計者にとっては作りが楽だけど、大容量で使う利用者が 対処しないといけないという安易なつくりなんですよね。 その代り少量データを対象とした統計手法研究などには 楽なんで。 Rにはなんといってもtidyverseがあるからなあ dplyrとかpurrrとか簡単で強力だから解析やるならR まあ適材適所で、ベイズならstanだし、機械学習や文字ならPythonが優勢、事務処理ならExcelVBAってとこ >>171 統計って回帰分析から急に難解になる それにじっと耐える生徒さんは偉い R初心者の私にはとっても役立つ講義だ 何かいいエディタないかな Rに最初からついてるエディタ使ってるけど不便 >>184 RStudioは2,3回インストールしてみたけど良く分からないし重いしですぐ消してしまってたんだよね 慣れると便利なの? RStudioのエディタは入力補完機能が凄い便利。パッケージの関数だけでなく自分で作った変数や関数も補完してくれる。 関数のヘルプもキー一つで表示できるし、ノートブック機能を使えばエディタ内で実行結果を表示してくれる。 その他、パッケージ管理とかプロジェクト管理とか便利だと思う。 具体的な話を聞くと便利そうな気がしてくる やっぱりRStudioがベストなのかもね しばらく使ってみようかな、ありがとう 補完機能なしでプログラミングするのって、武器を持たずに福岡の街を歩くようなものだよ。 初心者には Rコマンダーは便利だが 基本が理解するのが課題 R version 3.4.3 でR consoleを起動した直後に乱数を発生させると set.seedで指定していないのに前回起動時と同じ乱数が出てしまいます。たとえば runif(5)だと [1] 0.3108217 0.4556531 0.2194121 0.6496887 0.2677915 が毎回出てきます。 そういのものなのでしょうか? >>192 R x64 3.4.1(windows10)ではそれを再現できなかった そういうものではないと思う >>193 ありがとうございます。そうなんですか。自分の環境はwin7sp1 64bit で Rは5年前にインストールしその後何度か上書きインストールしたまに使っててたんですが 何らかの理由でRの中の設定が変更されたか破壊されてるんですかね 原因がわかったらまた書きます。 >>192 R3.4.3@Win10だけどそういう状態にはならないねえ。RStudioとかVSCodeとかでも何度かやってみたけど全て違う値。 事前にhtmlwidgetとか使ってない? >>192 です 試しに別のフォルダに新規インストールして何度か起動を繰り返して確認してみたら 正常に毎回違う乱数がでました。 なのでもう新規インストールで実行することにしました。 でも、どこに問題があるのか特定しないとまた起きそうなので少しづつパッケージを削除しつつ確認中・・・ >>195 ありがとうございます。 htmlwidgetsって何?って感じなんですが入ってました。 いつ何のために入れたのかもわからないのですが(笑) >>192 です。 パッケージの問題ではなくて、パッケージの作者の方々疑ってごめんなさい どうも>>195 さんの指摘どおり事前になにかが読み込まれてること気づきました 起動直後にカーソルの↑を押すとなぜか 昔コンソールで打ち込んだコマンドが未だに出てきちゃうんです たとえば、 > source('~/R/test.r') .Rhistoryがおかしくなっちゃったのかな?まあ、原因が分かってなにより? >>192 です。 原因がわかりました。 [Previously saved workspace restored]と出てることに今気づきまして 以前マイドキュメントに保存していたワークスペースが自動で読み込まれるためでした。 でもなぜか、"test.RData"の名前で保存したものは自動で読み込まれず ".RData"のように拡張子のみのファイル名が存在するとそれが読み込まれるようです。 新規インストールのRではがあっても読み込まれないのですが、 試しにC:\Program Files\R\R-3.4.3\bin\x64のフォルダ内にコピペしたら起きるようなので 変な名前やいけない場所には保存しては駄目ってことですね ありがとうございました。 .Rhistoryでset.seed()使ってるってことなのか >>200 うちのマシーンの環境がめちゃくちゃの可能性もありますが、 新規インストール直後のコンソールで保存したワークスペースをそれで読み込んでも 同じ乱数が再現されるので 保存した時のその後の乱数生成パターンも復元されると解釈しちゃいました。 はぇ〜 そういや最近出たブルーバックスのRの本が売れてるみたいだね >>201 環境の自動復旧は無効にしておいた方がいいよ。ググれば無効にする方法が見かるハズ… >>202 あの本、ブログか何かで宣伝してたよ。センセーショナルな見出しだったな。 この記事だ ttps://lab-on.jp/article/23 その本を買ってPCで真似してます Rコマンダーのお陰で敷居は下がりました Rの壁は 出力データの読み方、意味するところの理解 質のよいデータはあるので頑張る 共立出版のRで学ぶデータサイエンスシリーズの15巻とか16巻は ずっと欠けてるけどいつ出るのですか? >>204 jupyter notebook と pandas 最強だな エクセルソルバーのGRGみたいな非線形最適化制約有りをやりたいんだけど、Rでもできる? >>205 です 読み終えた、とても勉強になりました 困難さはRではなく統計学でしたが これは実践するしかないですね 気になったとこは頭に入れておきます P112の3次元散布図ですが Microsoft R Open version3.4(64bit)では 作図が表示されませんでした Rをかじっても やはり、考え方がわかってないとだめだね ということで、 「入門はじめての分散分析と多重比較」を読んでる Rだけとりあえず使えても、統計学がわかるわけじゃないからねw しかし慣れるとデータ扱うだけでも便利だなぁ。 でもちょっと不安なのが、変数が隠せないから自分で書いた変数とパッケージ内の関数が 偶然だぶったりしないもんなのかな・・・ってところ。 どうなん? ダブっても問題無いんじゃない? パッケージの変数を上書きしてしまう理由として考えられるのは 1. 知らずに上書きしてしまう 2. 知っててあえて上書きする だけど、どちらの場合も元の変数を使わないだけの気がする 安心するには実際にパッケージを壊してみるのも手かも 理由は異なるけど変数スコープの問題でwithを使ってる my = list (`%>%` = purrr::`%>%`, json = jsonlite::toJSON, id = function (x, ...) x); my = with (my, { id (1 : 10) %>% json %>% print; # my環境の変数を使う rocal = 123; rocal %>% print; # rocalは外から見えない my $ environment = function (...) 'hello world'; # 勧められないけど無問題 my $ try = tryCatch; # あえて上書き my; }); with (my, { try (typeof (rocal), error = function (ex) 'rocal is local') %>% print; environment (try) %>% print; # hello world base::environment (try) %>% print; # my$tryの環境はbase base::environment (try) = base::environment (); # my$tryの環境をmyに変更 tryCatch (try (stop ('try'), error = function (ex) 'catch') %>% print , error = id) %>% print; # my$tryが壊れた my$environmentが悪さ base::try (stop ('try')) %>% print; # base::tryは無傷 base::environment (try) = parent.env (base::environment ()); # my$tryの修理 try (stop ('try'), error = function (ex) 'catch') %>% print; # my$tryの復活 }); てな感じ 多くのサンプルでwithにデータフレームを代入してるけどリストでもオケー 小さな計算をいくつもする時にお手軽だと思う 今から学ぶならJuliaにせい Rよりずっと速い、書きやすい 僕の知り合いの知り合いができたパソコン一台でお金持ちになれるやり方 役に立つかもしれません グーグルで検索するといいかも『ネットで稼ぐ方法 モニアレフヌノ』 6O1Q6 >>136 あるよ ``` {python} def joking (a, b): while True: a, b = b, a + b; yield b; ``` ``` {r} with (new.env (), { obj = reticulate::py $ joking (0, 1); sapply (1 : 10, function (...) reticulate::iter_next (obj)); }); ``` もっとマシというより楽しい解答はこれかな 読んでないけど * [Python-like generators in R ・ GitHub](https://gist.github.com/klmr/d10623a0b4c7e1e9a6523eebee4913d1 ) で、最終的な答えは「ない」だと思う 技術的な問題というより文化的なところに理由がある気がする pythonにパイプ演算子がないのと同じ理由 実用的なyieldを実装しようとすると文法の修正が必要になると思う 言語仕様の複雑化という投資が必要になるので それなりの需要が見込めないと導入には踏み切れないかもね フィボナッチ数列ではyieldの有り難みは見えないけど 自動微分とは比較にならないほど汎用性の高い有り難い機能だと思う だけど自動微分と似ているところがあって どんな複雑なニューラルネットも泣けば自前で微分できるし どんな複雑な関数のyieldも泣けば自前でイテレーターを作れる fib = lambda n: int(((1+sqrt(5))/2)**n/sqrt(5) + .5) >>217 初めてききました。 調べたら6年前につくられてるようですが 書籍はほとんどないようですね。 なんで普及してないのでしょうか。 >>217 初めてききました。 調べたら6年前につくられてるようですが 書籍はほとんどないようですね。 なんで普及してないのでしょうか。 Juliaのこと? 本屋では単行本は1冊しか見かけないな。 あと、「データサイエンティスト養成読本 R活用編」というムック本の中に若干の記事があるくらいかな。 他にある? たぶん書籍が少ないのと、蓄積されたノウハウとか他人の作った関数とか少ないから まだみんな食いつかないんじゃないかな。 自分もなんか良さげではあると思いつつRで済ましてるしw 盛り上がっていたので貼っとく 読んでないけど Faster R with FastR | Hacker News https://news.ycombinator.com/item?id=18193557 FastRを使う人は少ないかもしれないけど 中のおしゃべりは楽しめるかも 個人的にはRotaさんの名前がここで出てきたことに驚き こういうところがHNの面白さかも Deriving the Normal Distribution | Hacker News https://news.ycombinator.com/item?id=18261892 正規分布からの連想で コメントにあるインタラクティブなデモに感心 A tutorial on Principal Component Analysis | Hacker News https://news.ycombinator.com/item?id=18256048 ブログの最後に書いてあるけど 現状はインタラクティブを実現するには 鬼プログラミングが必要なのかもしれない インタラクティブからの連想で Distill ― Latest articles about machine learning https://distill.pub/ ggmap使えなくなってた Googleにクレカ登録が必須なのか ggmapってなんだか知らないけどこの辺の話かな? Change in Google Geocoding API billing Issue #227 dkahle/ggmap GitHub https://github.com/dkahle/ggmap/issues/227 お買い上げありがとう御座います もう5年経ったか 機械学習データマイニングω流行ってるのに1/4スレも消費していない 次の3/4スレは15年じゃまだ余るだろう そしてリアルの理学部数学科の雰囲気に非常によく似ている。 Rのマニアルを見たら関数のパラメータがやから多くてびびった Rは難しいね =と<-が違うとか、forループ使うなとか Rは簡単やろ?forは禁止されてなく使っても構わない。 けどベクトル演算に持ち込んだほうが圧倒的に速い処理があるってだけの話。 プログラミングの基礎がないから、データフレームをnestしてmutateからのmap_dblという教科書的な処理すら少したつと忘れてしまうわ コンパイラー以外では余り語られないかもしれないパーサーの話 The new pqR parser, and R’s “else” problem https://www.r-bloggers.com/the-new-pqr-parser-and-rs-else-problem/ elseはC言語もしくはもっと古い言語の理論的バグと言ってもいいかも Dangling else - Wikipedia https://en.wikipedia.org/wiki/Dangling_else yacc/bisonの問題にエラーメッセージが理解不能なことがある Yacc is dead https://arxiv.org/abs/1010.5023 yaccそのものというよりbisonとの間で情報の欠落が起きるらしい Yacc is Not Dead (2010) | Hacker News https://news.ycombinator.com/item?id=8782218 パーサーを自動生成するのは人間がプログラミングするには複雑過ぎるからだけど 複雑 in 理解不能 out あれ?何処かで見たような おまけ:現在のニューラルネットの興隆は彼との対話が起点らしい Geoff Hinton Facts http://yann.lecun.com/ex/fun/ >232 斜め読みできる分量ではないけどオーサムらしい Free online book: Geocomputation with R, a book on geographic data analysis, visualization and modeling : rstats https://www.reddit.com/r/rstats/comments/a5tmht/free_online_book_geocomputation_with_r_a_book_on/ Geocomputation with R - the afterword | R-bloggers https://www.r-bloggers.com/geocomputation-with-r-the-afterword/ 各トピック毎に使えるソフトウェアの紹介があるので乗り換えの参考になるかも 中に書いてあるようにギットハブにもリソースがある GitHub - Robinlovelace/geocompr: Open source book: Geocomputation with R https://github.com/Robinlovelace/geocompr#geocomputation-with-r Geocomputation with R https://geocompr.github.io/ >>247 ありがとう こんな良さげな資料が公開されてたとは でもggmapのように地図(日本語)に図形やテキストをプロットした画像を出力する手法は載ってないみたいで残念 こういう話ではないと思うけど ``` {r} leaflet::leaflet () %>% leaflet::addTiles () %>% leaflet::addMarkers ( lng = 174.768, lat = -36.852, popup = 'The birthplace of R' ) %>% leaflet::addLabelOnlyMarkers ( lng = 174.768 + 0.01, lat = -36.852, label = '此処は何処ですか?' , labelOptions = leaflet::labelOptions ( noHide = T, textOnly = T, opacity = 1, style = list ( 'font-size' = '2ex', 'color' = 'red' ) ) ) %>% leaflet::addRectangles ( lng1 = 174.768 - 0.02, lat1 = -36.852 - 0.02 , lng2 = 174.768 + 0.02, lat2 = -36.852 + 0.02 ) ``` leafletではさすがに希望には合わないかな 仕事で自分のクレカ使うのも嫌だし、こういう有料化は困る タイトル勝ち - envのハッシュテーブルとしての使い方 Hash Me If You Can | R-bloggers https://www.r-bloggers.com/hash-me-if-you-can/ 分散の占有率は考えたことがなかった Principal Component Analysis (PCA) 101, using R : rstats https://www.reddit.com/r/rstats/comments/akytig/principal_component_analysis_pca_101_using_r/ PCAネタでOjaの学習則 - ご本人による解説 Oja learning rule - Scholarpedia http://www.scholarpedia.org/article/Oja_learning_rule スティーフェル多様体上の勾配 The Geometry of Algorithms with Orthogonality Constraints https://arxiv.org/abs/physics/9806030 動機と出発点は異なるが同一の微分方程式を導き出してる 対称行列の大きい方の固有ベクトルを求める問題は次のようにも書ける given M: (n,n)対称行列, to be found X: (n,k)行列 argmax <X, M X> subject to <X, X> = 1, where <A, B> := tr (A^T, B) subject to ...の部分がスティーフェル多様体の定義になっている 特に、k=1の場合は(n-1)次元単位球面になる 身近にある多様体の例になっていると思う 間違え subject to X^T X = (k,k)単位行列 ``` {Rcpp} #include <Rcpp.h> // 「rcpp 参照渡し」で検索すると出てくる話題 // [[Rcpp::plugins(cpp14)]] // [[Rcpp::export]] SEXP unsafe_negate (SEXP out) { // コピーなしを確実にするために面倒だがSEXP switch (TYPEOF (out)) { // 思いつく残りのキーワード: Rtools on Windows, case REALSXP: { // knitr::all_rcpp_labels, RCPP_MODULE, Rcpp/dispatch.h. Rcpp::NumericVector a (out); // RCPP_MODULE = boost::python std::transform (a.begin (), a.end (), a.begin (), std::negate <double> ()); return out; // https://wiki.python.org/moin/boost.python/HowTo } break; default: { // R 3.5以上で動くかわからない throw std::runtime_error ("unsupported type"); } break; // https://purrple.cat/blog/2018/10/14/altrep-and-cpp/ } // デビアン系だけかもしれないけど、Rstudioの環境下ではC++が超お手軽 } // [[Rcpp::export]] Rcpp::NumericVector safe_negate (Rcpp::NumericVector inn) { Rcpp::NumericVector out = Rcpp::no_init (inn.size ()); std::transform (inn.begin (), inn.end (), out.begin (), std::negate <double> ()); return out; } ``` 残りを貼り忘れた ``` {r} a = 1.0; b = safe_negate (a); cat (a, '->', b); # 1 -> - 1 a = 1.0; b = unsafe_negate (a); cat (a, '->', b); # - 1 -> - 1 ``` 統計使ってオプションの自動売買したいんですけど これ使えばできるようになりますかね? もう30年近く前からいっぱい使われてるよ。 それより統計知識と業務知識が先にないとあかんよ。 最近ハマったRの文法 ``` {r} tryCatch ({ f = function (x) { x + 1; } + 1; cat (f (1)); }, error = function (ex) { cat (ex); }); ``` 落とし穴というより他のプログラミング言語からの固定観念に縛られていた jsと対比してみる <pre id = 'dump'></pre> ``` {js} $ (window).on ('load', function () { try { const f = function (x) { x + 1; } + 1; $ ('#dump').html (f (1)); } catch (ex) { $ ('#dump').html (ex); } }); ``` 仕事で急遽Rを勉強している者です。 rvestを用いたスクレイピングについて解る方 いらっしゃいましたらご教授願います。 @下記ページを参考に rvest html_sessionで「次へ」のリンクを辿ってURLを抜き出し http://estrellita.hatenablog.com/entry/2015/11/11/084310 AそのURLをリストに追加して Bread_html をかけようとしてるんですが Cno applicable method for 'xml_find_all' applied to an object of class "list" とエラーが出てしまいます。 無知で申し訳ないのですが、原因と解決法わかりますでしょうか? わかる人が来るまでのつなぎ - 30年以内に来ればラッキーだと思うけど * ['rvest::html' is deprecated, but rvest::read_html doesn't exist. Issue #191 tidyverse/rvest GitHub](https://github.com/tidyverse/rvest/issues/191 ) * [rvest package | R Documentation](https://www.rdocumentation.org/packages/rvest/versions/0.3.2 ) ``` {r} with (list (`%>%` = purrr::`%>%`, size = length, null = NULL), { home = 'https://stackoverflow.com' ; depth = 0; done = list (); todo = list ('/questions/28863775/scraping-linked-html-webpages-by-looping-the-rvestfollow-link-function'); while (0 < size (todo) & depth < 2) { depth = depth + 1; done = c (done, todo); todo = purrr::reduce (.init = null, .x = todo, .f = function (out, path) { url = paste0 (home, path); text = xml2::read_html (url); nodes = rvest::html_nodes (text, css = '.related a.question-hyperlink'); purrr::reduce (.init = out, .x = nodes, function (out, node) { path = rvest::html_attr (node, 'href'); if (path %in% done) { cat ('what a small world:', path, '\n'); out; } else { cat ('i am going to stalk:', rvest::html_text (node), '\n'); c (out, path); } }); }); } }); ``` no applicable method for 'xml_find_all' applied to an object of class "list" あちこちの変数をデバッグすれば? まず、エラーの場所を特定するべき! スクレイピングには、Ruby, Nokogiri, Selenium WebDriver が良い なんで、こういうツールってDOMでのスクレイピングじゃないの? 覚えるの面倒なんだ。 バグ修正 if (path %in% done | path %in% out) { cat ('what a small world:', path, '\n'); out; } else { cat ('i am going to stalk:', rvest::html_text (node), '\n'); c (out, path); } 楽しいライブラリ * [Lego Mosaics Using R | Hacker News](https://news.ycombinator.com/item?id=19469142 ) ライブラリの階層としては [rgl](https://github.com/cran/rgl )/ [rayshader](https://github.com/tylermorganwall/rayshader )/ [brickr](https://github.com/ryantimpe/brickr )/ という感じかな? Rのライブラリというと変化球勝負というイメージを持っているのだけど rglは豪速球な気がする 誰でも頭が良くなる、プログラムが書けるようになる方法が発見される 81137 https://you-can-program.hatenablog.jp * [TF-IDF in a nutshell](https://www.reddit.com/r/LanguageTechnology/comments/bb5bcr/tfidf_in_a_nutshell/ ) [TF-IDF](https://en.wikipedia.org/wiki/Tf%E2%80%93idf ) が覚えられないので [PMI](https://en.wikipedia.org/wiki/Pointwise_mutual_information ) と関連付けてみる \newcommand{\nwd}[2]{\sharp\left\{{#1}\to{#2}\right\}} $X$を単語の有限集合、$Y$を文書の有限集合とする データを単語から文書への [二部グラフ](https://en.wikipedia.org/wiki/Bipartite_graph )として見て 辺の統計を考える | 記号 | 日本語の記号 | |:--|:----------------------------------| | $\nwd{x}{y}$ | 単語$x$の文書$y$での出現頻度 | | $\nwd{x}{*} := \sum_{y\in Y}\nwd{x}{y}$ | 単語$x$の全文書での出現頻度 | | $\nwd{*}{y} := \sum_{x\in X}\nwd{x}{y}$ | 文書$y$の長さ | | $\nwd{*}{*} := \sum_{x\in X,\; y\in Y}\nwd{x}{y}$ | 全文書の長さ | $$ \frac{P(x,y)}{P(x,*)P(*,y)} := \frac{\cfrac{\nwd{x}{y}}{\nwd{*}{*}}} {\cfrac{\nwd{x}{*}}{\nwd{*}{*}}\cfrac{\nwd{*}{y}}{\nwd{*}{*}}} = \underbrace{\frac{\nwd{x}{y}}{\nwd{*}{y}}}_{\approx\mathtt{TF}} \underbrace{\frac{\nwd{*}{*}}{\nwd{x}{*}}}_{\approx\mathtt{IDF}}. $$ * [March: "Top 40" New CRAN Packages](https://www.r-bloggers.com/march-2019-top-40-new-cran-packages/ ) [lenses](https://cran.r-project.org/web/packages/lenses/index.html ) というライブラリが入っている ゲット/セットのペアからスタートする [nlab](https://ncatlab.org/nlab/show/lens+%28in+computer+science%29 ) ではプットと書かれているがセットと同じ [イントロ](https://cfhammill.github.io/lenses/index.html ) も簡潔な説明だと思う 「てか、`iris $ Sepal.Length [3]`という書き方で何も困っていないし」 と思うかもしれないけどそれは正常な感覚だと思う [この例](https://ncatlab.org/nlab/show/Grothendieck+group ) はもっとビミョーな気分になるかもしれない 「ここで小麦粉を$1/3$カップ入れます」の$1/3$という書き方は [随伴](https://ncatlab.org/nlab/show/adjoint+functor )を表している 正の自然数からスタートした場合、Haskell風の書き方をして 直積に同値関係`Bunsu (a, b) == Bunsu (c, d) = a * d == b * c`を定義すると 単位射が`return a = Bunsu (a, 1)`で 積射が`join (Bunsu (Bunsu (a, b), Bunsu (c, d))) = Bunsu (a * d, b * c)`の モナドになる [この例](https://en.wikipedia.org/wiki/Simplicial_set ) はRの`list (list (0 : 0), list (0 : 1), ..., list (0 : n))`という リストがつくる圏からスタートする UMAPの土台に使われているらしい * [UMAP](https://cran.r-project.org/web/packages/umap/vignettes/umap.html ) * [How UMAP Works](https://umap-learn.readthedocs.io/en/latest/how_umap_works.html ) Rで正規表現を使って文字列から数値だけ抽出するにはどう書きますか? Ruby では、\d で数字を、+ で1文字以上、to_i で整数型に変換する p ary = "1a23bc04".scan( /\d+/ ).map( &:to_i ) #=> [1, 23, 4] ただし、これでは、負数を処理できない。 負数を処理するには、- を、? で、0か1文字 p ary = "-1a23bc-04".scan( /-?\d+/ ).map( &:to_i ) #=> [-1, 23, -4] 小数点や、e 表記は、もっと難しい。 これらは正規表現じゃなく、ライブラリを探すべき! どんな文字列か分からんけど単に数字を取り出すなら stringr::str_extract_allで[[:digits:]+]を指定すりゃいいのでは? パターンによってはreadr::parse_numberでいけるけど >>269 ありがとうございます。無事目的を果たすことが出来ました。 >>270 [[:digit:]+]だったね。いつも間違える。 最終的にどうやったか書いてくれると同じように困ってる人が助かるよ >>271 そうですね、では貼っておきます。 str <- "10.5万人" res <- as.numeric(stringr::str_extract_all(str,pattern = "[0-9.]+")) [The first web site](https://news.ycombinator.com/item?id=13356062 ) まで行くと古すぎだけど [Static Web - Back to the Roots](https://news.ycombinator.com/item?id=20006850 ) knitrでHTMLを作るときMathJaxをSVGに変換したくなることもある nodejsが動くようだと[MathJax-node](https://github.com/mathjax/MathJax-node )が簡単 Webネタで [Helping One Million Developers Exit Vim (2017)](https://news.ycombinator.com/item?id=19069526 ) Rを境に上位はWebのフロントエンドで下位は汎用プログラミング言語 つまり`R = (JQuery + Ruby) / 2` ウソだけどどんなウソにも真実がある気もする >>268 ま、Rubyで統計解析しなくてもいいんじゃね? 今のところ、信頼性保証できないでしょ PerlからRuby覚えてかなり気に入っていたけど、最近の機械学習ブームで 完全にPythonに負けた感があるな 今から始める人にRubyやれとは言えなくなった Python一択の時代 大学最後の2ヶ月で、元彼、好きだった人、女友達、男友達の合計500人の特性を全てデータ化して、関係の継続年数と好き度合いをアウトプットにして、Rでモデル組んだの。 それに結婚候補者の10人のデータいれて、1番結果良かったのが今の旦那で、付き合って結婚したの。本当に統計学は最強の学問。 [so]{#so}: [Showing existence of a diffeomorphism preserving volume forms](https://math.stackexchange.com/questions/2731058 ) 一次元だと積分表が使えるので絵を描いてみる [so-1]{#so-1} ``` {r} big_data = list (`%>%` = purrr::`%>%`, add = rlist::list.append , size = length, null = NULL, true = T, false = F, na = NA); big_data = with (c (big_data, xa = 0, sa = 0.1, db = 1, sbm = 0.3 , sbp = 0.2), { pa = function (x) { dnorm (x, mean = xa, sd = sa); }; pb = function (x) { 0.5 * (dnorm (x, mean = - db, sd = sbm) + dnorm (x , mean = db, sd = sbp)); }; Pa = function (x) { pnorm (x, mean = xa , sd = sa); }; Pb = function (x) { 0.5 * (pnorm (x, mean = - db, sd = sbm) + pnorm (x, mean = db, sd = sbp)); }; ra = function (n) { rnorm (n , mean = xa, sd = sa); }; rb = function (n) { sample (c (rnorm (n , mean = - db, sd = sbm) , rnorm (n, mean = db, sd = sbp)), n); }; x = 3 * seq (- 1, 1, len = 1e+3); ya = pa (x); yb = pb (x); plot (range (x) , range (ya, yb), type = 'n'); lines (x, ya, col = 'blue'); lines (x, yb , col = 'red'); add (big_data, pa = pa, pb = pb, ra = ra, rb = rb, Pa = Pa , Pb = Pb); }); no_plot = function (text = 'space') { plot (c (0, 1), c (0 , 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n'); text ( x = 0.5, y = 0.5, text);}; ``` 青の分布を赤の分布に連続的に変形させる [so-2]{#so-2} ``` {r} big_data = with (c (big_data, nx = 2e+3, nt = 1e+3, alim = 2, blim = 4), { vec = function (x, t) { (Pa (x) - Pb (x)) / ((1 - t) * pa (x) + t * pb (x)); }; t = seq (0, 1, len = nt); dt = diff (t) %>% mean; x = ra (nx); x = x0 = x [which (abs (x) <= alim)]; x = purrr::map (t, function (t) { x <<- x + dt * vec (x, t); }) %>% do.call (what = rbind); draw = function (t, x, nbin = 20 , vector = null) { xlim = range (x, na.rm = true); plot ( xlim, range (t), type = 'n', main = 'trajectories', xlab = 'x', ylab = 't'); purrr::map ( sample (1 : ncol (x), 50), function (j) { x = x [, j]; if ( all ( is.finite (x))) { lines (x, t); } else { points (x, t); }; }); if ( is.function (vector)) { no_plot (); }; doit = function (ind) { t = t [ind]; x = x [ind, ]; if (is.null (nbin)) { nbin = 'scott'; }; msg = sprintf ( 'at the time %.1e', t); msg = c (msg, sprintf ('%d / %d', sum ( xlim [1] <= abs (x) & abs (x) <= xlim [2], na.rm = true), size (x))); hist (x, breaks = nbin, freq = false, xlim = xlim, main = msg); x = seq ( xlim [1], xlim [2], len = 1e+3); # lines (x, pa (x), col = 'blue'); lines (x, pb (x), col = 'red'); if (is.function (vector)) { x = seq (min ( xlim), max (xlim), len = 1e+5); plot (x, vector (t, x) %>% abs (), type = 'l', log = 'y', main = msg); }; }; for (i in c (1, size (t) / 2, size ( t))) { doit (i); } }; x = rbind (x0, x [- nrow (x), ]); x [which (blim < abs (x))] = na; fix_bin = function (x, lim) { max (sum (abs (x [nrow (x) , ]) <= lim , na.rm = true) / 50, 10) %>% as.integer (); }; draw (t, x , nbin = fix_bin (x, blim), vector = vec); add (big_data, x0 = x0, t = t , dt = dt , nt = nt, nx = nx, vec = vec, draw = draw, fix_bin = fix_bin , alim = alim, blim = blim);}); ``` 逆方向に変形させてみる [so-3]{#so-3} ``` {r} with (c (big_data), { wec = function (x, t) { - vec (x, 1 - t); }; x = rb (nx); x = x0 = x [which (abs (x) <= blim)]; x = purrr::map (t, function (t) { x <<- x + dt * wec (x, t); }) %>% do.call (what = rbind); x = rbind (x0, x [- nrow (x), ]); x [which (blim < abs (x))] = na; draw (t, x, nbin = fix_bin (x, blim), vector = wec); }); ``` おしまい [so-2](#so-2)のバグ:誤: `vector (t, x)` 正: `vector (x, t)` でも誤の方が素直 - 数式上は`vector :: Time -> (Position -> Vector)` $$\newcommand{\calM}{\mathcal{M}}\DeclareMathOperator*{\arginf}{arg\,inf}$$ $\calM$を多様体とし、その上の [経験分布](https://en.wikipedia.org/wiki/Empirical_distribution_function ) $p:\calM\to(0,1)$が与えられた時、モデル$q:\Theta\to(\calM\to(0,1))$を [KL距離](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence ) $D$を使って、$\theta_*:=\arginf_{\theta\in\Theta}D(p,q_\theta)$と フィッティングするというのが典型的なパターンだと思う。[so](#so)では、 初期分布$\alpha:=q_0$と推定した分布$\beta:=q_{\theta_*}$を繋ぐ座標変換の 集合$\phi_t:\calM\to\calM$を具体的に作っている。写像$\phi_t$で$\alpha$が $q_t$になったとすると、<em>粒子数が保存すべし</em>という要請は次のように 書かれる。$$\int_{x\in\phi_t(D)} q_t(x) = \int_{x\in D} \alpha(x) \quad\text{for all}\quad D\subseteq\clM.$$$\phi_t(x)$が$x$について微分できる ことを仮定すれば、この式は$\phi_t^*q_t=\alpha$という形にまとめられ、 <em>$q_t$を$\phi_t$で [プルバック](https://en.wikipedia.org/wiki/Pullback_ (differential_geometry)) すると$\alpha$になる</em>と読む。$\phi_t(x)$が$t$について微分できることと、 $x$について可逆なことを仮定すると、次の式が得られる。$$(\partial_t +d\iota_{X_t})q_t = 0\quad\text{with}\quad\phi_t^*q_t = \alpha\eqtag{adv}$$ 導出の経緯から、この式は頻出問題になっていて、 [移流](https://en.wikipedia.org/wiki/Advection )という名前がついているが、 あまりに頻出過ぎて、多くの分野で名無しになっていると思う。問題は、 <em>$q_1=\beta$となる$\eqref{eq:adv}$の解があるか?</em>ということになる。 [so](#so)では、Moserのトリックを使って、そのような解があることを示している。 もう少し続ける。 Moserのトリックの副産物として、$\eqref{adv}$は [ガウスの法則](https://en.wikipedia.org/wiki/Gauss%27s_law )になる。 $\calM$が一次元の場合、電場$\iota_{X_t}q_t$が [CDF](https://en.wikipedia.org/wiki/Cumulative_distribution_function )に、 クーロンポテンシャルがCDFの積分に対応する。CDFのような単調増加関数を [シグモイド関数](https://en.wikipedia.org/wiki/Sigmoid_function )で近似する ことは自然なことだろう。同じことだが、 [LogSumExp](https://en.wikipedia.org/wiki/LogSumExp )でCDFの積分を近似する ことも自然なことだろう。さらに、コードのように、初期分布と推定した分布が 局在化している場合は、局在箇所にシグモイドを対応させることで、効率的に 良い近似が得られる。 コードでバグったのはベクトル場$X_t$をプロットしているところだが、$X_t$の プロットは、動作確認にはなるが、変化が激し過ぎて挙動の理解には向いていない。 $X_t$の計算では、性質上、ゼロ割に近い状況が発生することが避けらない。コード では、そのケアに追われている。一方、Moserのトリックのおかげで、電場は時刻に 依存しない。シミュレーションを通して成り立つ大雑把な挙動の把握には電場の方が 向いている。 神経にできることは森にもお隣さんにもベイズにも安藤モアにもできる、多分。 逆も真なり。これに<em>ガウスにもできる</em>が加わっただけかもしれないが、 今まで気が付かなかった。 [ありがとう](https://www.youtube.com/watch?v=Fd3_uMxG608 )、 [いきものがかり](https://arxiv.org/abs/1806.07366 ) おしまい [ggwaterfall](https://www.reddit.com/r/rstats/comments/c3rf6m/ggwaterfall_a_package_for_drawing_density_and/ ) を[DFT行列](https://en.wikipedia.org/wiki/DFT_matrix )をネタに使ってみる [参考](https://en.wikipedia.org/wiki/Chebotarev_theorem_on_roots_of_unity ) ``` {r} with (c (big_data, n = 17, N = 1e+3), { out = outer (1 : n, 1 : n, function (i, j) { ij = (i * j) %% n; k = pracma::gcd (ij, n); angle = 2 * pi * (ij %/% k) / (n %/% k); complex (real = cos (angle), imaginary = sin (angle)); }); print (eigen (Conj (out) %*% out, only = true) $ value / n); out = purrr::map ((n - 2) : 2, function (k) { purrr::map (1 : N, function (.) { out = out [sample.int (n, k), sample.int (n, k)] %>% eigen (only = true); # prod (out $ values); prod (out $ values) / sqrt (k); }) %>% unlist (); }); `%.%` = function (f, g) function (x) f (g (x)); purrr::map (out, log %.% Mod) %>% ggwaterfall::waterfall_ft (); # purrr::map (out, log %.% Mod) %>% ggwaterfall::waterfall_density (); }); ``` DFT行列のチェックのバグ修正+testthat+lenses ``` {r} with (c (big_data, n = 17, N = 1e+3, get = lenses::view, set = lenses::set), { `%.%` = function (f, g) function (x) f (g (x)); values_l = lenses::index ('values'); square = function (a) (t %.% Conj) (a) %*% a; out = outer (1 : n, 1 : n, function (i, j) { ij = (i * j) %% n; k = pracma::gcd (ij, n); angle = 2 * pi * (ij %/% k) / (n %/% k); complex (re = cos (angle), im = sin (angle)); }); testthat::test_that ('dft', { # testthat::expect_equal (square (out), n * pracma::eye (n)); testthat::expect_equal ((Re %.% square) (out), n * pracma::eye (n)); testthat::expect_equal ((Im %.% square) (out), pracma::zeros (n)); }); k = (n - 2) : 2; out = purrr::map (k, function (k) { purrr::map (1 : N, function (.) { out = out [sample.int (n, k), sample.int (n, k)]; out = eigen (out, only = true) %>% get (values_l); (sum %.% log) (out) - 1 / 2 * log (k); }) %>% get (lenses::unlist_l); }) %>% set (lenses::names_l, k); testthat::test_that ('det', purrr::map (out, function (out) { testthat::expect_true ((all %.% is.finite) (out)); })); purrr::map (out, Re) %>% ggwaterfall::waterfall_ft (show.labels = true); }); ``` バグ修正 誤:(sum %.% log) (out) - 1 / 2 * log (k); 正:(sum %.% log) (out) - k / 2 * log (k); 正規化の方法を修正する必要があるが DFT行列を[直交行列](https://en.wikipedia.org/wiki/Orthogonal_matrix ) に変えても似た挙動をする 関係ははっきりしないが [この話](https://www.r-bloggers.com/two-interesting-facts-about-high-dimensional-random-projections/ ) を思い出した 新たにバグを生み出した可能性が高いがレンズの話 可逆な行列$M$でパラメトライズされたゲットを$\mathtt{get}_M(x):=Mx$とすると レンズ則を満たすセットが$\mathtt{set}_M(x,y)=M^{-1}y$と一意に定まる [レンズの可逆性](https://www.twanvl.nl/blog/haskell/isomorphism-lenses ) の例になっていると思う レンズ則をチェックしてみる [setup]{#setup} ``` {r} big_data = list ( `%>%` = purrr::`%>%` , add = rlist::list.append , size = length , null = NULL , true = T , false = F , na = NA); big_data = with (c (big_data, get = lenses::view, set = lenses::set, why = T), { plain = function (M) { if (why) { # avoid lazy evaluation if (is.matrix (M) != true | nrow (M) != ncol (M)) { stop ('parameter must be an invertible matrix.'); } } lenses::lens ( view = function (state) { as.vector (M %*% state); }, set = function (state, value) { as.vector (solve (M, value)); } ); }; value_l = lenses::index ('value'); make_state = function (value) { list (hello = 'world') %>% set (value_l, value); }; make_lens = function (M) { lenses::`%.%` (value_l, plain (M)); }; add (big_data, make_state = make_state, make_lens = make_lens, get = get, set = set); }); ``` ``` {r, cache = F} with (c (big_data, nx = 3, nt = 1e+3), { normalize = function (x) x / sqrt (sum (x * x)); testthat::test_that ('lens', purrr::map (1 : nt, function (.) { a = pracma::randortho (nx) %>% make_lens (); s = rnorm (nx) %>% normalize () %>% make_state (); x = rnorm (nx) %>% normalize (); y = rnorm (nx) %>% normalize (); testthat::expect_equal (s %>% set (a, s %>% get (a)), s); testthat::expect_equal (s %>% set (a, x) %>% get (a), x); testthat::expect_equal (s %>% set (a, x) %>% set (a, y), s %>% set (a, y)); })); }); ``` [setup](#setup)で`why`を`F`にするとテストが通らない どこにバグがあるかわからないが 遅延評価のチェインで露呈したバグになっている おしまい 上のバグの再現 バージョンによるかもしれないし仕様なのかもしれない ``` {r} library (purrr); f_plain <- function (x) { function (y) { x + y; }; }; f_guard <- function (x) { x <- x; function (y) { x + y; }; }; tryCatch ({ x <- 1 %>% f_plain (); cat (x (2), 'success piped-plain\n'); }, error = function (msg) { cat ('failed piped-plain:\n'); print (msg); }); tryCatch ({ x <- 1 %>% f_guard (); cat (x (2), 'success piped-guard\n'); }, error = function (msg) { cat ('failed piped-guard:\n'); print (msg); }); ``` RstudioからSSHでポートフォアードさせてDBに繋ぎたいのですがどのように記述すれば良いのでしょうか 君が本気で情報を得たいのであれば 今やるべきことは煽ることではない 足りない情報を出すことだ ## Configuration over Configuration $$ \require{TeX/extpfeil}\require{TeX/AMScd} \newcommand{\eqtag}[1]{\tag{#1}\label{#1}} \newcommand{\isa}[1]{\mathinner{\left[\!\left[{#1}\right]\!\right]}} \DeclareMathOperator{\bbR}{\mathbb{R}} \DeclareMathOperator{\ecdf}{\mathtt{ecdf}} \DeclareMathOperator{\ebdf}{\mathtt{ebdf}} \DeclareMathOperator{\relu}{\mathtt{relu}} \DeclareMathOperator{\lure}{\mathtt{lure}} $$ ``` {r} big_data = list ( size = length , add = rlist::list.append , test = testthat::test_that , lty_none = 0 , lty_line = 1 , lty_dash = 2 , lty_dot = 3); prelu = function (...) pmax (..., 0); plure = function (...) pmin (..., 0); ``` 続く ## Convex on Rails 多分、[村人の定理](https://en.wikipedia.org/wiki/Mathematical_folklore ) 同値な関数達 ``` {r} ebdf_0 = function (xi) { n = length (xi); qi = cumsum (sort (xi)); function (x) { purrr::reduce (.init = 0, .x = 1 : n, .f = function (out, j) { pmax (out, j * x - qi [j]); }) / n; }; }; ebdf_1 = function (xi) { n = length (xi); function (x) { purrr::reduce (.init = 0, .x = xi, .f = function (out, xi) { out + prelu (x - xi); }) / n; }; }; ebdf_huge = function (xi) { n = length (xi); function (x) { purrr::reduce (.init = 0, .x = 1 : n, .f = function (out, j) { purrr::reduce (.init = out, .x = combn (1 : n, j, simplify = F), function (out, js) { pmax (out, j * x - sum (xi [js])); }); }) / n; }; }; ebdf_2 = function (xi) { n = length (xi); xi = sort (xi); qi = cumsum (xi); function (x) { purrr::reduce (.init = n * x - qi [n], .x = n : 1, .f = function (out, j) { out - plure ((out + qi [j]) / j - xi [j]); }) / n; }; }; dog_data = with (big_data, { equal = testthat::expect_equal; xi = c (- 1, 0, 2, 3); test ('ebdf', { doit = function (xi) { bdf_0 = ebdf_0 (xi); bdf_1 = ebdf_1 (xi); bdf_2 = ebdf_2 (xi); bdf_huge = ebdf_huge (xi); x = seq (min (xi) - 1, max (xi) + 1, len = 1e+3); equal (bdf_0 (x), bdf_1 (x)); equal (bdf_0 (x), bdf_2 (x)); if (size (xi) < 5) { equal (bdf_0 (x), bdf_huge (x)); } }; doit (xi); n = 10; doit (c (rnorm (n, - 1, 1), rnorm (n, 1, 2))); }); add (big_data, xi = xi, equal = equal); }); ``` 続く Rの関数`ecdf`は次のように定義されている。 $$ \begin{split} \ecdf(x) &:= \ecdf_\xi(x) := \frac{1}{n} \sum_{j=1}^n \isa{x \ge \xi_j}, \\ \isa{\operatorname{expr}} &:= \begin{cases} 1, & \text{ iff } \operatorname{expr} = \mathtt{true}, \\ 0, & \text{ otherwise}. \end{cases} \end{split} \eqtag{eq:ecdf} $$ ここで、$\xi_1<\cdots<\xi_n\in\bbR$を観測された値とする。 これを次のように積分したものを$\ebdf$と書くことにする。 $$ \ebdf(x) := \ebdf_\xi(x) := \int_{y=-\infty}^x \ecdf(x) = \frac{1}{n} \sum_{j=1}^n \relu(x - \xi_j). \eqtag{eq:ebdf} $$ `ebdf_1`はこの式を、`ebdf_0, ebdf_huge, ebdf_2`は式変形したものを実装している。 ``` {r} with (dog_data, { plot (ecdf (xi)); x = seq (min (xi) - 1, max (xi) + 1, len = 1e+3); bdf = ebdf_2 (xi); plot (x, bdf (x), type = 'l', main = 'ebdf'); points (xi, bdf (xi)); lines (x, prelu (x - mean (xi)), lty = lty_dash); }); ``` コードを見ると、`ebdf_2`は、 [ReLU](https://en.wikipedia.org/wiki/Rectifier_ (neural_networks))を [活性化関数](https://en.wikipedia.org/wiki/Activation_function )とする [ResNet](https://en.wikipedia.org/wiki/Residual_neural_network )に なっていることがわかる。さらに、`ebdf_2`を微分すると、ReLUとシグモイドが 混在したResNetになる。 ``` {r} ecdf_2 = function (xi) { n = length (xi); xi = sort (xi); qi = cumsum (xi); function (x) { purrr::reduce (.init = cbind (n * x - qi [n], n), .x = n : 1, .f = function (out, j) { f = out [, 1]; d_f = out [, 2]; g = (f + qi [j]) / j - xi [j]; cbind (f - plure (g), d_f - (g < 0)); }) / n; }; }; ``` 続く ``` {r} with (dog_data, { test ('resnet-cdf', { doit = function (xi) { cdf = ecdf (xi); cdf_2 = ecdf_2 (xi); x = seq (min (xi) - 1, max (xi) + 1, len = 1e+3); equal (cdf (x), cdf_2 (x) [, 2]); }; doit (xi); n = 10; doit (c (rnorm (n, - 1, 1), rnorm (n, 1, 2))); }); }); ``` 何が嬉しいのかはさておき、$\eqref{eq:ecdf}$がResNetで書けたことになる。 おしまい ## Zigzag on Rails [ブログ](http://blog.mirkoklukas.com/finite-sample-expressivity/ )の内容を 実装してみる。 ``` {r} zigzag_on_rails = function (xi, yi) { js = order (xi); xi = xi [js]; yi = yi [js]; y0 = yi [1]; ai = diff (c (0, diff (yi) / diff (xi), 0)); f = function (x) { purrr::reduce (.init = y0, .x = seq_along (xi), .f = function (out, j) { out + ai [j] * prelu (x - xi [j]); }); }; big_data $ add (big_data, xi = xi, yi = yi, y0 = y0, ai = ai, zigzag = f); }; dog_data = with (dog_data, { yi = c (2, 0, - 1, 1); rails = zigzag_on_rails (xi, yi); x = seq (min (xi) - 1, max (xi) + 1, len = 1e+3); plot (x, rails $ zigzag (x), type = 'l'); points (xi, yi); add (dog_data, yi = yi); }); ``` 観測値からの外延の仕方がブログとは異なると思うが、こんな感じじゃないかと 思う。[区分線形関数](https://en.wikipedia.org/wiki/Piecewise_linear_function ) を式にすると、自然とReLUの和が現れる。折れ線の変化`ai`を正と負の成分に 分けると、それぞれが凸 on Railsに乗る。 続く ``` {r} convex_on_rails = function (xi, ai) { js = order (xi); xi = xi [js]; ai = ai [js]; bi = cumsum (ai); ci = cumsum (ai * xi); f = function (j, x) bi [j] * x - ci [j]; d_f = function (j, x) ai [j] * (x - xi [j]); i_f = function (j, x) (x + ci [j]) / bi [j]; convex = function (x) { js = seq_along (xi); purrr::reduce (.init = 0, .x = js, .f = function (out, j) { pmax (out, f (j, x)); }); }; resnet = function (x) { js = rev (seq_along (xi)); purrr::reduce (.init = f (js [1], x), .x = js, .f = function (out, j) { out - plure (d_f (j, i_f (j, out))); }); }; list (convex = convex, resnet = resnet); }; resnet_on_rails = function (xi, yi) { out = zigzag_on_rails (xi, yi); with (out, { xi = out $ xi; ai = out $ ai; js = which (ai > 0); fp = convex_on_rails (xi [js], ai [js]); js = which (ai < 0); fm = convex_on_rails (xi [js], - ai [js]); convex = function (x) y0 + fp $ convex (x) - fm $ convex (x); resnet = function (x) y0 + fp $ resnet (x) - fm $ resnet (x); add (out, convex = convex, resnet = resnet); }); }; ``` 続く ``` {r} with (dog_data, { rails = resnet_on_rails (xi, yi); x = seq (min (xi) - 1, max (xi) + 1, len = 1e+3); plot (x, rails $ resnet (x), type = 'l'); points (xi, yi); doit = function (xi, yi) { rails = resnet_on_rails (xi, yi); x = seq (min (xi) - 1, max (xi) + 1, len = 1e+3); equal (rails $ zigzag (xi), rails $ convex (xi)); equal (rails $ zigzag (xi), rails $ resnet (xi)); }; doit (xi, yi); n = 10; doit (c (rnorm (n, 0, 1), rnorm (n, 2, 3)), c (rnorm (n, 0, 2), rnorm (n, 0, 1))); }); ``` 区分線形関数は二本のResNet with ReLUの [アフィン写像](https://en.wikipedia.org/wiki/Affine_transformation )で 書けたことになる。 おしまい と思ったけどバグを見つけた 誤: equal (rails $ zigzag (xi), rails $ convex (xi)); equal (rails $ zigzag (xi), rails $ resnet (xi)); 正: equal (rails $ zigzag (x), rails $ convex (x)); equal (rails $ zigzag (x), rails $ resnet (x)); いいと思う だけどここは過疎を通り越して廃墟なので返事は期待できないと思うよ >>307 タイムスタンプ >>1 が20年くらい前かと思ったが 意外と新しくて驚いた 素晴らしい [dataAnim](https://github.com/chrk623/dataAnim ) ``` {r} with (list (), { no_thanks = function (...) data.frame (..., stringsAsFactors = F); j = 1 : 3; a = no_thanks (name = LETTERS [j], a = j, b = letters [j]); j = j + 1; b = no_thanks (name = LETTERS [j], a = j, c = letters [j]); dataAnim::join_anim (join_type = "left", speed = 1 , x = a, y = b, by = 'name', show_msg = T); }); ``` [JavaScript const in R](https://colinfay.me/js-const-r/ ) やってみる ``` {js} $ (window).on ('load', function () { const out = new Map (); out.hello = 'world'; $ ('#const-here').html (JSON.stringify (out)); }); ``` <pre id='const-here'></pre> ``` {r} with (list (), { your_name = function (x) deparse (substitute (x)); out = new.env (); lockBinding (your_name (out), environment ()); tryCatch ({ out $ hello = 'world'; }, error = function (ex) { assign ('hello', toString (ex), env = out); }); jsonlite::toJSON (as.list (out), auto = T); }); ``` ``` {Rcpp} #include <Rcpp.h> // [[Rcpp::export]] void kossori (Rcpp::List out, Rcpp::String key, Rcpp::RObject val) { out [key] = val; } ``` ``` {r} with (list (), { your_name = function (x) deparse (substitute (x)); out = list (hello = ''); lockBinding (your_name (out), environment ()); tryCatch ({ out $ hello = 'world'; }, error = function (ex) { kossori (out, 'hello', toString (ex)); }); jsonlite::toJSON (out, auto = T); }); ``` 仕様上は、`env`に対する`$`でのセッターは、たとえ`lockBinding`が 施されていても無問題と思うが、実際は、弾くようにしている。 何か理由があるのかな? [ブログ](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系関数は、カラム名に日本語があると正しく処理されないバグがないか? 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 ) こうした言葉を思いつくのは才能なんだろうね。 Rだと簡単にハイゼンバグの例を作れた。 ``` {r a_10850, dependson = ""} uncertainly = with (new.env (), { plus = \(x) \(y) x + y x = 1L plus_1 = plus (x) x = 1000L plus_1 (2) |> print () }); ``` アレレ?ということでデバッグ文を入れてみる。 ``` {r a_29042, dependson = ""} principle = with (new.env (), { plus_debug = \(x) { cat ("Schr?dinger says ", x, "\n", sep = "") \(y) x + y } x = 1L plus_1 = plus_debug (x) x = 1000L plus_1 (2) |> print () }); ``` デバッグ文の有無で挙動が変わる。`uncertainly`の挙動は、バグではなく、 言語仕様だと思った方が良いと思う。 `uncertainly`の挙動を防ぐには関数`force`を使えば良い。 * [How to not fall into R's 'lazy evaluation trap](https://stackoverflow.com/questions/29084193 ) タイトルに"遅延評価"が入っているが、遅延評価自体の問題ではなく、 Rの実装に原因がある。Haskellも遅延評価だが、`uncertainly`のような 挙動が起きれば、Haskellはその存在意義を失う。 ``` ghci import Data.Functor.Identity (Identity (..)) :{ uncertainly :: (Monad f) => f Int uncertainly = do x <- pure 1 plus_1 <- pure (x +) x <- pure 1000 pure $ plus_1 2 :} runIdentity uncertainly flip ($) () uncertainly maybe 0 id uncertainly ``` <pre> 3 3 3 </pre> 言語仕様レベルでのトレードオフかもしれない。 * 局所最適化と大域最悪化 * 整合性とコンパイル時間 * 等など magrittrの例を少し改変 # fns first <- function(x){ message("first") invisible(x) } second <- function(x){ message("second") invisible(x) } # lazy NULL %>% first() %>% second() # eager NULL %!>% first() %!>% second() ビックリパイプは知らなかった。"抜かりなし"だね。ただし、ビックリパイプは ハイゼンバグとは関係ないかな。 ハイゼンバグから離れて、宇宙大戦争について書いてみたい。 `lapply`等のループ系の関数は、最終的にC関数`R_forceAndCall`を呼び出す。 * [apply.c](https://github.com/wch/r-source/blob/trunk/src/main/apply.c ) Rは、"フォースと共にあらんことを"などと呑気なことを言っている場合ではなく、 フォースと共にあらねば死んでしまう。 前回のコード`uncertainly`は変数`x`の使い回しが敗因だが、ループ系の関数 では変数を使い回すしかない。そのために、ループ系の関数ではフォースが必須 になっている。98パーセントぐらいの使用例では、フォース抜きでも動作する (当社調べ)が、残りの2パーセントでコケる。僅か2パーセントでも、 プログラミングをギャンブルにしないためには、穴を塞いでおく必要がある。 上のコード`uncertainly`の挙動を"プロミスの罠"と書くことにする。 ここでの"プロミス"は、JSの"プロミス"ではなくて、Rでの"変数"の 実装方法を指す。 * [6 Functions | Advanced R](https://adv-r.hadley.nz/functions.html ) この記事には次の一節がある。 > You cannot manipulate promises with R code. Promises are like a **quantum state:**: ... 多分、ここでの"量子状態"はハイゼンバグと同じ現象を指しているんだと思う。 観測すると、状態が変化してしまう。 関数`do.call`自体はプロミスの罠と関係しないと思う。 ``` {r a_28697, dependson = ""} do.call (`+`, list (1L, 2L)) ``` Pythonだと次のコードに対応する。 ``` {python a_10786, dependson = ""} (lambda x, y: x + y) (* range (1, 3)) ``` 殆どのプログラミング言語で、関数の引数リストは [一級市民](https://en.wikipedia.org/wiki/First-class_citizen ) でないように思う。一級市民でない代わりに、一級市民の配列からの変換が 用意されている。Pythonでは`*`という関数がビルトインで用意されている。 Rでは引数リストを直接作れるかもしれないが、`do.call`で配列から 変換するのが一般的だと思う。 自分の知る限り、プロミスの罠にハマるのは次のパターンに限られる。 ~~~ {.r} x = "hello" g = f (x) x = "world" g () ~~~ "関数を返す関数"`f`に、"変数"`x`を代入した時にプロミスの罠が可能性が 出てくる。Rの評価戦略はHaskellと同じ [コールバイニード](https://en.wikipedia.org/wiki/Evaluation_strategy ) に分類されている。上の例では、次の場合にプロミスの罠が発生する。 1. 関数`f (x)`の中で引数`x`がニードにならず、 1. 返り値の関数`g`に渡される。 次の例はプロミスの罠が発生する。 ``` {r a_20914, dependson = ""} rude = with (new.env (), { const = \(x) \(...) x x = "hello" g = const (identity (x)) x = "world" g () |> print () }); ``` 関数への代入はニードでないので(by definition)、関数適用の連鎖 `const (identity (x))`の中にはニードがない。そのために、最終的に`g ()`が コールされるまで引数`x`は評価されない。コード`g ()`がコールされた時点で 引数`x`の定義を探しに行くので、プロミスの罠にハマる。 コールバイニードを実現するためには、コンパイラーが内部的に次のように 書き換える必要があると思う。 ``` {r a_1451, dependson = ""} polite = with (new.env (), { const = \(x) \(...) x x_1 = "hello" g = const (identity (x_1)) x_2 = "world" g () |> print () }); ``` [静的単一代入](https://en.wikipedia.org/wiki/Static_single-assignment_form ) と呼ばれる操作と同じだと思う。この書き換えはループ処理では必須になる。 コード`rude`はアカン奴として切り捨てることができるかもしれないが、 ループ処理では変数の多重定義が避けられない。Rでは、何らかの理由で、 この書き換えが省略されているために、プロミスの罠が発生するのかな?と 思っている。 コード`rude`からコード`polite`への書き換えを上とは異なる形で行ってみる。 Rではできないので、Pythonを使う。 ``` {python a_14213, dependson = ""} ##| cache: false # def none (): id = lambda a: a const = lambda x: lambda * y: x def lhs (): x = "hello" g = const (id (x)) x = "world" return g () def rhs (): return (lambda x: (lambda g: (lambda x: g ()) ("world")) (const (id (x)))) ("hello") assert lhs () == rhs () none (); ``` 関数`lhs`から関数`rhs`への書き換え規則は次のようになる。[rule]{#rule} ~~~ {.python} y = f (x) z = g (x, y) rest (x, y, z) == (lambda y: z = g (x, y) rest (x, y, z) ) (f (x)) == (lambda y: (lambda z: rest (x, y, z)) (g (x, y))) (f (x)) ~~~ この書き換え規則を適用した結果、関数`rhs`では、静的単一代入への変換は、 [ラムダ計算](https://en.wikipedia.org/wiki/Lambda_calculus )のアルファ変換 に置き換わる。この意味では、静的単一代入とアルファ変換は同じことになる。 実際のコードでは、分岐、ループ、副作用が入ってくるので、こんな単純は 話では済まないと思うが、極度に単純化すると、静的単一代入とアルファ変換は、 単に、方言の違いということになる。 Rの場合は、書き換え規則[rule](#rule)は等価な変換にならない可能性がある。 特に、プロミスの罠が現れた場合は、挙動が変わると思う。 余談だが、書き換え規則[rule](#rule)をモナドに拡張したものがHaskellの [do-記法](https://en.wikibooks.org/wiki/Haskell/do_notation )になっている。 ~~~ {.python} y <- f (x) z <- g (x, y) rest (x, y, z) == kleisli-extension (lambda y: z <- g (x, y) rest (x, y, z) ) (f (x)) ~~~ モナドが恒等関数の場合に[rule](#rule)に一致するので、do-記法を [rule](#rule)の拡張と言って差し支えないと思う。 おしまい [トランスフォーマー](https://en.wikipedia.org/wiki/Transformer_ (machine_learning_model)) についての記事: オーサーヒアーかな? * [Coinductive guide to inductive transformer heads](https://news.ycombinator.com/item?id=34970877 ) 中を読んでいないが、アテンション機構を [ホップ代数](https://en.wikipedia.org/wiki/Hopf_algebra )で説明、もしくは、 ホップ代数を使ってアテンション機構の類似物を導出するという話のようだ。 この手の救世主的な話はがっかりすることが多いので、過大な期待は禁物かも。 g = do.call(const, list(x)) ではどうだろう RのS5クラス定義(setRefClass)ってsave、loadでシリアライズ、デシリアライズ出来ないのか デシリアライズされたメンバ関数にアクセスできない ## Rにおけるオブジェクト指向への取り組み 趣旨は新規フレームワークR7の紹介だが、既存のフレームワークの俯瞰図 としても使いやすいと思う。 * [What is R7? A New OOP System for R](https://www.jumpingrivers.com/blog/r7-oop-object-oriented-programming-r/ ) ## オブジェクト指向とは? 次の記事を取り上げてみる。 * [Object Oriented Programming Features of Rust](https://doc.rust-lang.org/book/ch17-00-oop.html ) 内容については賛否両論だと思う。しかし、箇条書きできるということが、 オブジェクト指向が単一の概念でないことを示しているように思う。 ## S3はオブジェクト指向か? 比較のために、S3のキーワード`default`を使わないで書く。 ``` {r a_3770, dependson = ""} none = with (new.env (), { to_string = \(x, ...) UseMethod ("to_string", x) as_string = \(x, ...) UseMethod ("as_string", x) to_string.Base = \(x, ...) as_string (x, ...) as_string.Base = \(x, ...) to_string (x, ...) to_string.Derived = \(x, ...) paste0 ("hello ", class (x)) as_string.Derived = as_string.Base structure (identity, class = "Derived") |> as_string () |> print () }); ``` `Base`クラスは、関数`to_string`か関数`as_string`の"どちらか1つを実装せえ"と言っている。 似た内容をPythonで書いてみる。 ``` {python a_8560, dependson = ""} class Base: def to_string (self): return self.as_string () def as_string (self): return self.to_string () class Derived (Base): def to_string (self): return "hello " + type (self).__name__ Derived ().as_string () ``` Haskellで書いてみる。 ``` haskell class Base_ a where to_string :: a -> String to_string = as_string as_string :: a -> String as_string = to_string instance Base_ () where to_string :: () -> String to_string = ("hello " ++) . show as_string () ``` <pre> "hello ()" </pre> ワンセットの関数群をインターフェースと書く。 1. Pythonのクラスでは、 * 構造体の定義とインターフェースの定義と実装が同時に行われる。 1. HaskellのクラスやRのS3では、 * 構造体とインターフェースは別々に定義される。 * インターフェースの実装は既存の"構造体"に対して行われる。 HaskellやRの"構造体"には、通常の構造体に加えて、関数も含まれる。 構造体とインターフェースが別々に定義される点で、RのS3はHaskellのクラスに 近いように思う。 おしまい join_by()が便利で感動した 気象データのマージがすごくシンプルになった 任意の地点・期間における測定値についての風向風速などの気象の影響をみるため近傍アメダスの気象データをマージする際に、便利になったと実感しました 具体的にはマージの際にbetween(時刻,)の条件が追加できるようになったのが大きいです 処理が重いかコードが長いかだった上のマージ処理が、join_by(地点名, between(time, start_time, end_time))で済むようになりました [非等価ジョイン](https://en.wikipedia.org/wiki/Relational_algebra )は 最近になって関数`dplyr::xxx_join`に導入されたらしい。 * Rで非等価結合 (2) * [dplyr 1.1.0: Joins](https://www.tidyverse.org/blog/2023/01/dplyr-1-1-0-joins/ ) [ラッキー](https://www.youtube.com/watch?v=xBJMr1v5Zuw )。 関数`vctrs::vec_locate_matches`には、以前から限定された形の非等価ジョイン のオプションがあったので、実用的な形にするのに時間がかかったのかもしれない。 リスト内包表記が使いたいのでPythonで書く。 ``` {python a_32372, dependson = ""} #| cache: false # def slow_matches (pred, zipper): #{ def go (a, b): #{ a = pandas.DataFrame ({"" : a}).groupby ("").groups b = pandas.DataFrame ({"" : b}).groupby ("").groups key = ((j, k) for j in a.keys () for k in b.keys () if pred (j, k)) val = (zipper (a [j], b [k]) for (j, k) in key) val = zip (* itertools.chain (* val)) return (numpy.fromiter (val, dtype = numpy.int32) for val in val) #} return go #} def none (): #{ a = numpy.random.choice (8, 2 * 8, replace = True) b = a + 4 def go (pred, zipper): #{ (j, k) = slow_matches (pred, zipper) (a, b) out = pandas.DataFrame ({"a" : a [j], "b" : b [k]}) print (out) #} go (lambda a, b: a == b, itertools.product) go (lambda a, b: a == b, lambda a, b: itertools.islice (zip (a, b), 1)) go (lambda a, b: abs (a - b) <= 1, itertools.product) go (lambda a, b: abs (a - b) <= 1, lambda a, b: itertools.islice (zip (a, b), 1)) #} none (); ``` `dplyr::xxx_join`の数々のオプションは理解できていないが、`NA`絡み以外は、 コールバック`pred`と`zipper`でかなりカバーしてるんじゃないかと思う。 ジョインの操作自体は明朗会計だが、実行時性能がボッタクリ価格に転嫁される というプログラミング社会が抱える構造的な問題が反映されている。 おしまい 一応Rスレだしコード一式を載せるならRにしてほしいな 正直、書き方を思い出すためだけに、話題にかこつけて、わざと異なる プログラミング言語で書くこともある。しかし、今回は純粋に、処理の道筋を 単刀直入な形で表せるプログラミング言語を選択した。 [線形回帰](https://en.wikipedia.org/wiki/Linear_regression )を例にとる。 Rで線形回帰はとても書きやすい。しかし、線型回帰をRのイディオムとして 覚えるよりも、線形代数として覚えた方が長い目で見た時に利益になると思う。 同じように、等価ジョインは [プルバック](https://en.wikipedia.org/wiki/Pullback_ (category_theory)#Sets) に対応する。等価ジョインを、Rのイディオムとして覚えるよりも、 単純な`for`ループの形で覚えた方が長い目で見た時に利益になると思う。 リスト内包記法は`for`ループを簡潔に書くための記法になっている。 プルバックの話の続きをアップした。 * [JSFiddle](https://jsfiddle.net/ytkhdpcq/show ) プルバックのカリー化の話になっている。プルバックの話はこれでお終い。 洗練された話は次の記事にある。 * [locally cartesian closed category in nLab](https://ncatlab.org/nlab/show/locally+cartesian+closed+category ) 表計算の立場で見ると、"何が問題なのかわからない"問題かもしれない。 ある規則を満たす表でサマリーするとプルバックのカリー化になる。 4.3.0入れてみた パイプがまたひとつ便利になっていい感じ Windowsでパスが長くても通るようになったのも嬉しい どなたか教えてくだされ 対照実験の医療統計をしたいと思って 傾向スコアマッチングにしようと思ってます 肺がん手術を実施した人たち200症例のうち10人だけ 手術後に合併症を起こしていました この場合 「合併症を起こしたグループ」・・・10人 「合併症を起こさなかったグループ」・・・190人になり 患者の背景で傾向スコアマッチングをキャリパー幅0.2でペア分けすると 各群8人ずつの比較になってしまいます。 これで優位差を出す統計するには数が少なすぎるのであまりにも乱暴な統計になりますか? それとも元々のサンプルが200人から選りすぐりのマッチングさせた8人ずつを選んでいるので8人ずつと少なくても問題ないですか? 共変量のバイアスがより少ないサンプルが得られますが、サイズが8人の二群比較なのは変わらないんじゃないかと思います >>493 つまりサイズが8人ずつの比較では 優位差を言うには弱いと思いますか? >>494 検定はサンプルサイズ込みの判定なので8人であっても二群に有意な差があったと言って差し支えありません 有意水準が甘ければ弱いとは感じますが 一方で、検定は標本のバイアスが除去できていることを保証してくれません スレッドの質をガクッと下げるが、Rの非標準評価で遊んでみた。 * [JSFiddle](https://jsfiddle.net/fo1qjt87/show ) 4.3.0のWindows版はdir()にバグあるね リストがまともに取得できない 4.2.3に戻した tibbleを利用してデータ探索しているんだけど、2つの変数で条件つけて散布図つくるにはどんな手順踏んだらいいのか… 変数1をプロットの色にして変数2をサイズにするとかじゃだめなの? それも一つの方法ですね。ありがとう ただ、目的以外の条件をプロットしたくないのです。 tibbleの構成としては、観測が4時点、変数がたとえば6個(うち1つは個体識別番号で、変数としてはA-Eの5種類あるとしましょう)のような感じで、 散布図のx軸に「観測時点1における全個体の変数A」 y軸に「観測時点4における全個体の変数E」 をプロットしたいのです。 個体識別番号をキューにして、異なる時点の異なる種類の変数を散布図にするとでもいいましょうか… その要件だと以下のような感じでいけると思います library(tidyverse) df <- tibble( time = rep(1:4, each = 5), id = rep(1:5, 4), a = rnorm(20), b = rnorm(20), c = rnorm(20), d = rnorm(20), e = rnorm(20) ) df |> pivot_wider(id_cols = id, names_from = time, values_from = c(a, e)) |> ggplot(aes(a_1, e_4)) + geom_point() ありがとうございます! それをヒントに試してみます。 >>507 要件の理解が違っていなければ散布図作成まで行きますよ time1でのaの値とtime4でのeの値をid毎にプロットしています >>508 失礼しました こちらのビューワの問題で、ブラウザでみたらスクリプトの最後まで表示されました💦 >>508 教えてもらったスクリプトでドンピシャでした、ありがとうございます。 どなたか統計素人の私に教えてください 医療でA群100名とB群100名をいろんな観点から 比較するとします 背景を揃えるために傾向スコア分析で キャリパー0.2で取って30ペア(合計60人)を抽出しました 例えば ①アミラーゼの値とかを統計にかけたいときは、その60人のアミラーゼのデータをまずは正規分布かどうかをShapiro-wilk検定で測って、0.05を上回っていたら「正規分布」と見なして、今度はこれら60人のアミラーゼのデータが等分散かどうかを確認するためにf検定を行なって等分散だったらt検定(Student's t-test)、違ったらWelch's t-testで優位差があるかどうかを調べる。 もしShapiro-wilk検定が0.05未満で非正規分布だったときはMann-whitney's U検定で優位差があるかどうかを調べる。 ↑Q1. この認識であってますか? Q2. この正規分布を計るのはペアを作った後のデータ(n=60)で、正規分布かどうか?等分散かどうか?を見ますか? それともマッチング前のデータ(n=200)で正規分布かどうか?等分散かどうか?を確認するべきですか? Q3. アミラーゼじゃなくて他の連続変数の項目(BMIだったり、血圧だったり、白血球数だったり)の優位差を調べる場合も 全て上の流れで一つ一つの項目ごとで正規分布か?等分散か?など確認していって適宜、該当する算出法を項目ごとに採用して優位差を測る必要がありますか? もし良ければ教えていただけるとありがたいです 素人だと思ううちは傾向スコアには手を出さないほうがいいような気がする >>512 上記の場合はどういうかんじでやるのが 最適解になるのかだけでも教えてもらえないでしょうか? writexlパッケージで保存すると時刻のタイムゾーンは反映されないのね ggraggedパッケージがなかなか便利 facet_grid()を詰めて並べられる R4.3.1にした dir()関連が直ってるようでひと安心 標準パイプの新機能もようやく使える ggplot2のカラーパレットのデフォルトが変わったりしました? 棒グラフの色の割り当てが変わった気がする… >>518 変わってないはず ただggplot2のデフォルトはグループ数に応じて関数で色を割り当ててるのでグループ数が異なると違うパレットのように感じるとは思います 株価データを取得したいがためにRを弄り始めたけど中々難しいね ググってコピペしたコードを一つ一つパッケージと関数の挙動を調べているけど やりたいこと全てができるようになるまでは遠そうだ >>520 データをどこから持ってくるかによって随分難易度は違う。 ファイルとしてダウンロードできるところもあるし、 スクレイピングが必要なところもある。場合によってはせれにうむとか。 企業が求めるエンジニア言語ランキング 2位「R」 https://www.itmedia.co.jp/news/spv/2306/20/news183.html 3位以下を大きく引き離しててすごいけど… Rユーザの少なさを反映しただけな気もする >求人数のトップ3は1位「C言語」(39万2559件)、2位「R言語」(30万3759件)、3位「Go」(8万3922件)だった。 >このランキングは、SalesNow DBに1月1日から5月31日までに載った114万5294件の有効求人から算出した。 5ヶ月で30万件もあったらしい Pythonが4万件弱なところを見るとデータ解析用言語を列記したわけじゃなく純粋にRの需要があったみたいだ 不思議 結局企業で欲しいのは統計処理だから、統計処理ならRだろうっていう短絡的な 求人が多かったってことじゃない? いくつか求人サイト見てみたけどRのみなんて求人はほぼ皆無でほとんどがPythonとの併記だった これ集計ミスかもね 同じ会社の別の期間のランキング見るとRは10位以内に入ってないし 業∞!!!! 字∞!!!!! 苗∞!!!!!! 米∞!!!!!!! 畑∞!!!!!!!! 灰∞!!!!!!!!! 糞∞!!!!!!!!!! 机∞!!!!!!!!!!! 結構何年か前にRは結構古臭いし、言語仕様もわかりづらいし衰退するのでは、と書いているのを 読んだことがあるけど、逆に今だと他の商用統計パッケージが廃れる勢いなのにRって立場が盤石な 気がする なんで? AIブームでPythonに全部もってかれると思ったけどそうでもなかったね 先発優位で統計屋や大学教授が昔からRユーザーってのもあるのかも? ユーザー視点からだとtidyverseの存在が大きいと思う 505のように地味に手間がかかりそうな処理が関数ひとつでひょいと処理できるとかすごい >>529 今でもC/C++が盤石なのと同じ理由な気がする とくに大学の研究室とかは、同じような機能のプログラムを別言語に コンバートする作業は雑用でしかないので誰も積極的にはやらない。 RでできるのならこのままRでいいじゃん、もしくはフォートランでもw みたいな感じ やったところでなにか新しい知見が得られるわけもなし。 昔の医学分野の統計処理はSASがデフォだったけど COVID-19で疫学の人(西浦先生とか)がRガンガン使ってるの見て疫学以外でもRが広がってる 学会行くとggplotで書いたグラフいっぱい見る ggplot2のデフォルトテーマは微妙なので変えているけどいつも悩む ggthemesをもとにそこから装飾することが多い ついにExcelでPython使えるようになるらしい マイクロソフトはRのベンダー買収してMRANとかやってたのに時代の流れかな Pythonでexcel弄るのは 15年以上前から出来てるのを今更 excel上のpythonで実行する方が遅い可能性まである PythonからExcelをいじることとExcelでPythonが動くことは全く違う pythonは何故かPCでJupyter notebookが不具合で使えなくてな 導入時点でコケるとその言語をやる気失くす…Rはその点、導入が簡単で良かった Excelも最近はTidydata志向で、パワークエリやパワーピボットで再現性のある加工や分析が できるようになってきているけど、使い始めはなかなかいいかな、と思ったがちょっと凝ると やっぱり複雑化してメンテできないVBAみたいになっちゃうなw ま、RやPythonだからメンテ性がいいというわけでもない気はするけどw Excelマクロ禁止してきた企業はPython機能も禁止しそう Rを社員が使うケースなんてのはもうレア度SSRなんだろうね Rで学ぶなんとか分析とかいう本いっぱい出ているし、結構個人でやる気のあるやつは覚えるんじゃね? R人口増えている印象。 昨日の統数研のワークショップ見て、大学や企業によってはライブラリのインストールに制限あったりでデフォルトで解析できるのがRが多言語に置き換わらない理由の一つなのかもなと思った 今後も「選択と集中」続けるどころか懸賞金型推進するようだから商用に金払うのも辛くなるとこも出てきそうね 数が増えてることと数が少ないこととは違う これ混同する人は結構多い 自分は今の部署にいる間だけは統計解析するけど異動したらもう解析なんてしないからRのモチベは低いな 人事の巡り合わせでいつかまた今の部署に戻るかもしれないけど 前回のおーまんちょーでは一部のコードで途中から行番号が表示されるようになっていた。 多分、JSFiddleのエディターで`tidy`というボタンを押してしまったからだと思う。 実は、`prettydoc::html_pretty`で行番号が表示されないのを直したいと思っていた。 しかし、どこから手を付けたら良いかわからなかった。それが、ハプニングのおかげで 手の入れ場所がわかった。多分、HTMLとしては不整合な状況になっているが、 意図した通りに表示されるようになった。 [ありがとう、おーまんちょー](https://jsfiddle.net/105nohbq/show ) group_by()はmap()と組み合わせて使えないんだね 一般的な意味では、そんなことないと思うよ。 状況によってはgroup_byとmapを組み合わせると上手く行くことがあると思う。 ぐるぐる図のくだりでグダグダなところがあったので、 [修正した](https://jsfiddle.net/6z7u3d09/show )。 アップした後にタイポを見つけたが、覆水盆に返らず。 多変数の時系列データをクラスタリングするpackageをご存知であれば教えていただけますか…? 知らない。 人がいれば、意味のない返事はしない方が良いかもしれないが、過疎だからね。 あなたは地球上で何日過ごしてますか? という問いに答えるRのスクリプト 干支と誕生日(月日)から今日が生まれて何日めか計算 eto=("ひつじ") tag="10-17" Eto=c("子(ね=ねずみ)","丑(うし)","寅(とら)","卯(う=うさぎ)","辰(たつ=龍(りゅう))","巳(み=へび)","午(うま)","未(ひつじ)","申(さる)","酉(とり=にわとり)","戌(いぬ)","亥(い=いのしし)") # Eto i=grep(eto,Eto) this_year=as.numeric(strsplit(as.character(Sys.Date()),'-')[[1]][1]) j=this_year%%12 - 3 age=12*(0:9) + j-i age=age[18<=age][1:5] age AD=paste0(this_year-age,"-",tag) days=\(x) as.numeric(Sys.Date() - as.Date(x)) rbind(age,days=sapply(AD,days)) あってますか? マジで聞いているんなら、テストツールの使い方を覚えた方が建設的だと思う。 個人的には、Haskellからの流れで [quickcheck](https://github.com/RevolutionAnalytics/quickcheck ) しか経験したことがないが、他にもあると思う。 コードが合っているかどうかわからないが、パット見問題なさそうに見える。 問題があったとしても、デバッグで十分対応できるように見える。 [リーナスの法則](https://en.wikipedia.org/wiki/Linus%27s_law ): あたなが十分な暇人ならばバグはいずれ塞がれる。 ChatGPTはなんていい奴なんだろう。 クソ野郎: I'd like to test the following code: ``` r 拝借 ``` * The code is written in the programming language R. * "子", "丑", ..., are the Chinese / Japanese zodiac called Eto in Japanese. This code have to calculate the age for given birth date in Eto, e.g. `eto = "ひつじ", tag = date = "10-17"`. Are there any advice to test this code? ChatGPT: Certainly! To test the provided R code, you can follow these steps: 1. Understand the Code: * Make sure you understand the logic of the code and how it calculates age based on the Eto and birthdate information. 1. Run the Code: * Copy and paste the code into an R script or R console. * Ensure that there are no syntax errors or issues when running the code. 1. Test Cases: * Define several test cases with different combinations of Eto, birthdate (tag), and current date. * For example, you can use different values for eto and tag, and verify the calculated age against your expectations. 1. Expected Results: * Manually calculate the expected age for each test case based on the provided Eto information and current date. * Compare the manually calculated age with the age calculated by the R code. 1. Edge Cases: * Test the code with edge cases, such as birthdates close to the boundary of age ranges, to ensure the code handles them correctly. 出禁になった。書けるかな? てなわけで、ちゃっとじーぴーてーの残りの部分は割愛する。 意地悪クソ野郎と親切じーぴーてーにも共通することろがあって、 最後の締めはやっぱりリーナスの法則になる。言い方が全然違うけど。 親切じーぴーてーが例に挙げたテストコードがテストの基本形になると思う。 くいっくちぇっくのようなテストツールでは、テストケースをランダムに発生させるが、 そうしたテストツールが使えない場合も多々ある。テストツールが使えない場合は、 親切じーぴーてーが示したようなベタな方法で行くしかない。 >>566 書き込めるかどうか判らぬままに 黙々とその低能感溢れる文章を書いている姿を想像すると泣けてくる 書き込めてよかったな 親切じーぴーてーとの会話を見返すと、関数を細分化して、(ナウ、干支、月日) という三つ組を引数とする関数を作るのもありかな?という気がしてきた。 このように関数を細分化すると、テストケースが作りやすくなる。 クラスタリングとか次元削減に関する教科書的な本ってありますか? かけるかな 統計の人でも機械学習の人でもないので、返事を書くのは心苦しいものがあるが、 過疎なのでしょうがない。統計処理と機械学習の区別がつかないぐらいの素人だよ。 ありきたりだけど: * [朱鷺の杜Wiki](http://ibisforest.org/index.php?FrontPage ) グーグルとアヒルで検索してみたけど、この検索については英語の方が情報が取りやすいと思った。 英語でも日本語でも本屋の広告が多くなるが、英語の方が俺が推薦図書はこれだ! みたいなのが多いような気がした。検索で見つけた本だけど非売品と商品の両方がある。 * [Modern Dimension Reduction](https://arxiv.org/abs/2103.06885 ) * [Modern Dimension Reduction](https://www.cambridge.org/core/elements/abs/modern-dimension-reduction/D1506D1D7D51CC4EF375E4633D54F377 ) このようなケースでは立ち読みしてから買うこともできる。時々、非売品と商品の 両方が用意されていることがあるので、注意した方が良いかもしれない。 殊、次元削減で言えば、線形代数みないなもので、ありとあらゆる場所で使われるので、 俺の次元削減にジャストフィットするものを探すのは難しいかもしれない。 ちなみに、専門的には、学校の教室で研究会を開くのが統計で、リゾート地で 国際会議を開くのが機械学習らしい。 read.cgi ver 07.5.1 2024/04/28 Walang Kapalit ★ | Donguri System Team 5ちゃんねる