統計解析R たぶんpart3くらい
>>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); }); ``` read.cgi ver 07.5.0 2024/04/24 Walang Kapalit ★ | Donguri System Team 5ちゃんねる