統計解析R たぶんpart3くらい
楽しいライブラリ * [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になるような read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる