【R言語】統計解析フリーソフトＲ第6章【GNU R】 [無断転載禁止]©2ch.net

**１３２人目の素数さん** · 2019/01/13(日) 19:54:34.41

>>649
知ってるかもだけどいちいちmappingは書かなくてもいい
scale_shapeも値が1:3ならいらない

ggplot(data, aes(shape=Conditions, col=Conditions))+
geom_line(alpha=0.6)+
geom_point(alpha=0.8)

**１３２人目の素数さん** · 2019/01/27(日) 04:45:13.06

習ってから思ったけど、Fortranとgnuplotでいいよね

**１３２人目の素数さん** · 2019/01/29(火) 23:19:36.86

普及度、パッケージ数、専門度、文献数の点でそれはない

**１３２人目の素数さん** · 2019/01/30(水) 01:27:02.50

ヒカキンの年収が10億超え!?明石家さんま・坂上忍も驚愕の総資産とは??
https://logtube.jp/variety/28439
【衝撃】ヒカキンの年収・月収を暴露！広告収入が15億円超え！？
https://nicotubers.com/yutuber/hikakin-nensyu-gessyu/
HIKAKIN(ヒカキン)の年収が14億円！？トップYouTuberになるまでの道のりは？
https://youtuberhyouron.com/hikakinnensyu/
ヒカキンの月収は１億円！読唇術でダウンタウンなうの坂上忍を検証！
https://mitarashi-highland.com/blog/fun/hikakin
なぜか観てしまう！！サバイバル系youtuberまとめ
http://tokyohitori.hatenablog.com/entry/2016/10/01/102830
あのPewDiePieがついに、初心YouTuber向けに「視聴回数」「チャンネル登録者数」を増やすコツを公開！
http://naototube.com/2017/08/14/for-new-youtubers/
27歳で年収8億円　女性ユーチューバー「リリー・シン」の生き方
https://headlines.yahoo.co.jp/article?a=20170802-00017174-forbes-bus_all
1年で何十億円も稼ぐ高収入ユーチューバー世界ランキングトップ10
https://gigazine.net/news/20151016-highest-paid-youtuber-2015/
おもちゃのレビューで年間12億円！今、話題のYouTuberは6歳の男の子
https://www.businessinsider.jp/post-108355
彼女はいかにして750万人のファンがいるYouTubeスターとなったのか？
https://www.businessinsider.jp/post-242
1億円稼ぐ9歳のYouTuberがすごすぎる……アメリカで話題のEvanTubeHD
https://weekly.ascii.jp/elem/000/000/305/305548/
世界で最も稼ぐユーチューバー、2連覇の首位は年収17億円
https://forbesjapan.com/articles/detail/14474
ヒカルの収入が日収80万、月収2400万、年収3億と判明www
https://matomenewsxx.com/hikaru-income-8181.html
はじめしゃちょーの年収は6億？2017年は30億突破か？
https://2xmlabs.com/archives/1873

**１３２人目の素数さん** · 2019/01/30(水) 11:36:05.44

glmer関数にgamma分布を適用したモデルをggplotで回帰曲線を引く方法がわかりません......。summary(model)$dispersionがNULLと返されてしまいます。

**１３２人目の素数さん** · 2019/01/30(水) 22:41:24.83

コードの例示がないからよく分からないけど、もともと戻り値が無いんじゃないの

**１３２人目の素数さん** · 2019/01/31(木) 17:27:54.13

?ggplotを見たかい

**１３２人目の素数さん** · 2019/02/02(土) 11:02:49.05

dplyr 0.8.0てもう来た？
2月1日とか聞いたような

**１３２人目の素数さん** · 2019/02/02(土) 21:51:39.00

>>656
可読性を保つために省略しないほうがいい
プログラミングの基本

**１３２人目の素数さん** · 2019/02/02(土) 22:39:28.49

逆じゃね。可読性上げるために省略する。
特にmappingなんかは書くまでもない。

**１３２人目の素数さん** · 2019/02/03(日) 00:29:26.31

継承を無視してコピペしまくるのが基本なわけない
マウント取りもいきすぎると滑稽だわ

**１３２人目の素数さん** · 2019/02/14(木) 03:33:35.52

>>666
お前のほうが滑稽
マウントってるのお前
イキってるだけで中身もないし
滑稽すぎる

**１３２人目の素数さん** · 2019/02/14(木) 07:38:25.36

データテーブルの1列目のある文字列（例えばstart）を検索してその行含めた上の行を全部削除するスマートな書き方ある？
df［-1:-grep("start",df[,1])］
今こんな感じだけどパイプで繋げない。
slice使えばいいのか、grep以外にいい関数があるのか？

**１３２人目の素数さん** · 2019/02/14(木) 09:16:23.85

>>668
startが必ず一つしかないとか、restartとかがあったときにそれを検出しても良いならそれでもよいが、普通は
df[-(1:which(df[, 1] == "start")[1]), ]
すると思う。
経験的にはgrepはできるだけ避けた方がいい。

**１３２人目の素数さん** · 2019/02/14(木) 11:48:21.71

お邪魔します。
スレ違いが明らかなんですが、ここで聞かせて下さい。

numpyスレ立てていいですか?

**１３２人目の素数さん** · 2019/02/14(木) 12:08:18.97

もちろん
とんでもなくスレ違い

さっさと立てるね

**１３２人目の素数さん** · 2019/02/14(木) 12:28:21.65

>>669
なるほど、ありがとうございます。
そしてカッコでくるんでからマイナス付ける書き方もあるのか。勉強になります。

**１３２人目の素数さん** · 2019/02/14(木) 12:45:38.37

>>667
可読性について勉強できてよかったな

**１３２人目の素数さん** · 2019/02/14(木) 17:44:07.85

データフレームのある列が全てNAの時、その列を削除するよい方法ある？
現状col<-apply(2,function(x)｛all(is.na(x))｝
で要らない列定義してからdf[,!col]としてるけどほんとはパイプの中に入れて処理したい。

**１３２人目の素数さん** · 2019/02/14(木) 23:27:05.62

>>674
自己解決しました。select_ifでいけそうです。

**１３２人目の素数さん** · 2019/02/16(土) 12:45:31.76

>>670
板違いだと思うので、やめておいた方が良いのでは

**１３２人目の素数さん** · 2019/02/19(火) 23:32:31.86

tidymodelsどう？あまり日本語の資料ないからなあ

**１３２人目の素数さん** · 2019/02/27(水) 15:36:28.78

RStudioをvi風のキーバインドにすると、
ノーマルモードのときに日本語入力してしまうとバグみたいになるんだが
あれどうにかならんの?

**１３２人目の素数さん** · 2019/02/28(木) 13:53:56.86

どなかた↑よろしく

**１３２人目の素数さん** · 2019/03/01(金) 22:53:08.03

RStudio と日本語入力ってほんと相性わるい

なんとかしてくれんかなぁ

**１３２人目の素数さん** · 2019/03/02(土) 11:44:51.96

windowsのpreview版(1.2.1303)RStudio使ってるけど、
IMEが無効になるバグが直ってる感じがする。
まだ十分使い込んでないからだけかも。

**１３２人目の素数さん** · 2019/03/03(日) 01:15:23.41

https://wired.jp/2019/01/18/get-wired-kevin-kelly-5-videos/
シロンボヒトモドキゴキブリニホンザルの自由は偽自由と詐欺広告人を殺す自由ヒトモドキゴキブリシロンボアメ公はニホンザルゴキブリと自殺せよ？

**１３２人目の素数さん** · 2019/03/04(月) 09:44:29.47

｢いきる｣とか最近ネットで使ってるやつ増えてきたがなんなんあれ?
あほな反抗期の中学生がつかってるイメージしかわかんのだが。

**１３２人目の素数さん** · 2019/03/04(月) 12:02:26.52

「生きる」と「熱る」のどっちだ？
いずれにしろ広辞苑に載ってるから調べればいい

**１３２人目の素数さん** · 2019/03/04(月) 12:42:02.72

【人類を２つに分ける】　世界教師マⅰトレーヤ到来
http://mevius.5ch.net/test/read.cgi/seiji/1538271437/l50

まもなく世界経済破綻、みなさんの財産は紙くずになります。

**１３２人目の素数さん** · 2019/03/04(月) 19:13:04.98

>>684
方言に対して広辞苑を持ち出すのは的外れ。
大阪弁辞典によると「意気る」または「粋る」なので、漢字も外れ

**１３２人目の素数さん** · 2019/03/04(月) 19:19:22.60

>>686
どっちも「熱る」の当て字やん

**saga** · 2019/03/05(火) 02:07:35.96

活きる

**１３２人目の素数さん** · 2019/03/05(火) 08:37:04.71

＞「生きる」と「熱る」

「熱る」って初めて聞いた
これは普通に使われているの？
方言？
生まれてこのかた、東京なんだけど聞いたこと見たことないんだが
単なる私の不勉強かな？

**１３２人目の素数さん** · 2019/03/05(火) 08:59:52.15

https://dictionary.goo.ne.jp/jn/10539/meaning/m0u/
いき・る【熱る】
１あつくなる。ほてる。むしむしする。
２激しく怒る。

https://kotobank.jp/word/%E7%86%B1%E3%83%BB%E7%86%85-2005439
いき・る【熱る】
① あつくなる。ほてる。むしむしする。
② 息づかいを荒くして怒る。相手と争おうとしていきまく。言いたてる。
③ 調子に乗って勢いこむ。元気づく。

**１３２人目の素数さん** · 2019/03/06(水) 20:09:12.61

https://www.youtube.com/channel/UCKyNitDuOcc7Gi1cjQlNCwA

ヒトモドキ奇形なりすましシロンボアメ公はイスラム国に首切られろヒトモドキフィフィラクダイスラム国の売女とともに自殺しろテロゴキブリ民族アメ公

**１３２人目の素数さん** · 2019/03/07(木) 21:17:38.50

豚肉屋の豚肉自民のキモオタ奴隷豚障害者犯罪者窃盗犯山田太郎と犯罪者キモオタは今すぐ民族レベルで自決自殺しろw
キモ豚山田太郎は自民入党のマイノリティ下僕化しキモオタヒトモドキ障害者の存在価値は0に達したんだよ表現戦士の性犯罪者害虫
このキモオタヒトモドキネトウヨをガス室で抹殺してえよな？w

**１３２人目の素数さん** · 2019/03/07(木) 22:45:56.15

腐れシロンボテロアメ公ヒトモドキは核で根絶やしになれヒトモドキニホンザルゴキブリの親玉障害者欠陥遺伝子白塵ゴキブリ

**１３２人目の素数さん** · 2019/03/10(日) 10:47:48.29

QAi2Acx7tz
奇形ネトウヨヒトモドキゴミ藪部落の顔気持ち悪いから自殺しろ糞犬hk

**１３２人目の素数さん** · 2019/03/10(日) 10:49:53.12

ww.sankei.com/column/ 180307/clm1803070005-a.html

ヒトモドキ産経便所ゴキブリは下着ドロの常習犯変態暴論ゴミ便所自殺しなさい今すぐ

**１３２人目の素数さん** · 2019/03/20(水) 12:05:11.55

欠損値のないデータで解析したモデルをstepAIC()で処理しようとしたところ、TRUE/FALSEが必要なところが欠損値ですとエラーが出たのですが、どう処理すればよいのでしょうか。MuMInのdredge()でも同様になってしまいました。

**１３２人目の素数さん** · 2019/03/28(木) 15:23:23.97

https://www3.nhk.or.jp/news/html/20190327/k10011863181000.html

インフルエンザの新しい治療薬「ゾフルーザ」を投与されたＡ香港型のインフルエンザ患者30人を調べたところ、70％余りに当たる22人から、この薬が効きにくい耐性ウイルスが検出されたことが国立感染症研究所の調査で分かりました。
調査件数は多くないものの、専門家は現在のような使用を続けると、耐性ウイルスが広がるおそれがあるとして使用基準を見直すべきだと指摘しています。

耐性化率が50%以上である確率は
pbeta(0.5,1+22,1+8,lower=F)
[1] 0.9946631でいいかな？

**１３２人目の素数さん** · 2019/03/28(木) 20:29:02.09

> (p=1-binom.test(22,30)$p.v)
[1] 0.9838752
> binom.test(22,30,conf=p)

Exact binomial test

data: 22 and 30
number of successes = 22, number of trials = 30, p-value = 0.01612
alternative hypothesis: true probability of success is not equal to 0.5
98.38752 percent confidence interval:
0.5000000 0.8994036
sample estimates:
probability of success
0.7333333

**１３２人目の素数さん** · 2019/03/29(金) 21:03:39.81

「統計」は「疑似科学」な

**１３２人目の素数さん** · 2019/04/06(土) 19:22:34.74

yahoo apiで距離取得するコード掲載したサイトないですかね
google料金高過ぎて鞍替えなんですが

**１３２人目の素数さん** · 2019/04/16(火) 20:17:02.49

>>680
使ってるQtが古いんや。1.2でだいぶましになった感じだけど入力途中の文字が残る妙な挙動が時々でる。

**１３２人目の素数さん** · 2019/04/16(火) 21:58:14.98

R書籍漁ったけどまともなのないな

**１３２人目の素数さん** · 2019/04/17(水) 07:23:57.02

>>702
書籍は出版される頃には情報が古くなってるから。自分はネットの情報とヘルプだけで十分だな。

**１３２人目の素数さん** · 2019/04/18(木) 19:43:25.70

名前が一文字のせいで検索する時に苦労する

**１３２人目の素数さん** · 2019/04/18(木) 22:13:27.80

つ seekR

**１３２人目の素数さん** · 2019/04/29(月) 14:15:34.67

RStudioって32ビット版って存在しないの？？

ダウンロードしようとしても32ビット版見当たらないし、Windows7,10用の64ビット版セットアッププログラムを起動しようすると32ビットへの選択とかはなくて、単にエラーになっちゃうし

誰か教えてください

**１３２人目の素数さん** · 2019/04/30(火) 14:31:58.91

httpstwitter.com/shotkr16

低脳中卒万引きヒトモドキネトウヨ猿ヒトモドキをさしころせ

**１３２人目の素数さん** · 2019/04/30(火) 17:29:05.11

>>706
古いのはここにあるけど今後のことをかんがえたらOSを64bit版に変えた方がいいよ。
https://support.rstudio.com/hc/en-us/articles/206569407-Older-Versions-of-RStudio

い · 2019/04/30(火) 20:32:55.24

ここでRをつかった仕事をされてる人いますか？
学生ですか？

仕事あれば教えてほしい

**１３２人目の素数さん** · 2019/05/01(水) 19:38:51.19

>>709
Rを仕事にしてはいないけど、仕事の一環でR使ってる
信頼性、入手しやすさ、扱いやすさのバランス的に自分の環境ではR一択なんだよな

**１３２人目の素数さん** · 2019/05/01(水) 23:10:05.28

>>710
ちなみにPythonはどう？

**710** · 2019/05/02(木) 08:41:31.17

>>711
AI目当てでこれから使おうと思って勉強中
今の職場で求められるのは主に統計用途や作図なんでRは手放せないかな

**１３２人目の素数さん** · 2019/05/02(木) 13:30:24.26

twitter.com/tukuhae

ゴキブリネトウヨヒトモドキ奇形売春婦肉便器なつこババア滅多刺しにして解体しろ

**１３２人目の素数さん** · 2019/05/02(木) 13:37:47.73

>>712
なるほど。
自分も機械学習やってるけど、R の方がやりやすいよ。ディープラーニングならPythonかもしれないけど

**１３２人目の素数さん** · 2019/05/02(木) 18:56:28.17

データ分析とか統計解析とかだとRかPythonかだね
機械学習よりは統計寄りだとRは仕事でもスタンダード、特にマニアックな統計解析だとRにしかパッケージがない
まあ使い慣れてるならRで機械学習やっても全然構わないけど

**１３２人目の素数さん** · 2019/05/03(金) 01:30:52.83

>>708
サンクス
まあメモリが2Gなので

**１３２人目の素数さん** · 2019/05/03(金) 01:33:11.23

>>715
少し前はRなんてだんだん廃れるみたいな雰囲気あったけど今はむしろ盛り返してデファクトな何時になってるからこの先10年は行きそうだなあ

機械学習特需でpythonの躍進はすごかったけど

**１３２人目の素数さん** · 2019/05/03(金) 18:06:51.34

2Gって四則演算もままならないだろ

**１３２人目の素数さん** · 2019/05/11(土) 16:40:31.81

Rstudio初心者です。
データファイル名を変えてフォルダには変更名で保存されてるのに、読み込むともとの名前のままなんです。どうすれば変更後の名前になるのでしょう？

**１３２人目の素数さん** · 2019/05/12(日) 20:48:36.39

>>708
サンクス
当面は32ビット版で頑張る

**１３２人目の素数さん** · 2019/05/17(金) 22:45:05.59

>>719
全く状況がわからん。
主語や目的語、必要な修飾語を省略せずに、分かるように説明しないと誰も助言できない。

**１３２人目の素数さん** · 2019/05/19(日) 05:50:12.72

いつのまにかどうでもいい国の資格試験の選択科目にまでなったしな

なんなんだろう

**１３２人目の素数さん** · 2019/05/19(日) 12:44:13.70

資格を馬鹿にする者は資格に泣く

**１３２人目の素数さん** · 2019/05/20(月) 08:40:12.19

資格にしがみついたらいつのまにか技術の進展から取り残された国がありましたとさw

**１３２人目の素数さん** · 2019/05/22(水) 10:42:04.74

>>719
それはsave, loadよ話か？

**１３２人目の素数さん** · 2019/07/11(木) 19:34:28.29

Pivot_longer & wider

**１３２人目の素数さん** · 2019/07/20(土) 11:06:35.86

0645
ふうL@Fu_L12345654321
学コン1傑いただきました！
とても嬉しいです！

ｈｔｔｐｓ://pbs.twimg.com/media/D-IuUuqVUAALnAB.jpg
ｈｔｔｐｓ://twitter.com/Fu_L12345654321/status/1144528199654633477
https://twitter.com/5chan_nel (5ch newer account)

**１３２人目の素数さん** · 2019/08/06(火) 14:27:52.67

> 臨時地震板
> 【備えあれば】防災用品非常食スレ121【憂いなし】
> https://mao.5ch.net/test/read.cgi/eq/1562388281/

**１３２人目の素数さん** · 2019/08/10(土) 11:28:08.57

ＲってMac版とWindows版で優劣ないの？
今から始めるならどっちがオススメ？

**１３２人目の素数さん** · 2019/08/10(土) 13:44:59.32

R自体に優劣はないが文字コード絡みの問題が多いWindows環境はおすゝめしない。

**１３２人目の素数さん** · 2019/08/10(土) 17:55:47.97

Macはフォントで苦労するけどな

**１３２人目の素数さん** · 2019/08/10(土) 21:06:01.89

じゃ、Linuxだ。

**１３２人目の素数さん** · 2019/08/10(土) 21:23:29.54

フォントで苦労ってどういう種類の？
文字化け？

**１３２人目の素数さん** · 2019/08/12(月) 04:58:06.03

>>729
R本体だとグラフィックディバイスに違いはあるけど、実質的に差はない。
でもRは本体だけで使うことは稀で、多種多様なGUIがあるから、
そこでOSの違いによる差が生じる。

**１３２人目の素数さん** · 2019/12/09(月) 08:29:17.97

面白い問題スレにあったのでシミュレーションしてみた。

# サイコロ

# 正６面体のサイコロがある．４面は青色、２面は赤色である．
# このサイコロを合計20回振るとき、最も起こりそうな順番はどれか？
# １．赤青赤赤赤
# ２．青赤青赤赤赤
# ３．青赤赤赤赤赤

sim <- function(){
a=sample(0:1,20, replace=TRUE, prob=c(4,2))
b=as.character(a)
c=paste(b,collapse="")
s1=paste(c(1,0,1,1,1),collapse="")
s2=paste(c(0,1,0,1,1,1),collapse="")
s3=paste(c(0,1,1,1,1,1),collapse="")
res=c(grepl(s1,c),grepl(s2,c), grepl(s3,c))
return(res)
}
k=1e6
re=replicate(k,sim())
mean(re[1,])
mean(re[2,])
mean(re[3,])

結果は、直感とおり、１が再頻

> mean(re[1,])
[1] 0.124672
> mean(re[2,])
[1] 0.080873
> mean(re[3,])
[1] 0.040564
>

grep使わない方法ってあるかな？

**１３２人目の素数さん** · 2019/12/10(火) 18:31:14.56

これって確率分布的には何になるんですっけ？

**１３２人目の素数さん** · 2019/12/11(水) 10:39:16.38

歪んだコインやんか

**１３２人目の素数さん** · 2019/12/12(木) 16:40:23.20

>>735
ワイの直感的解法

# １．赤青赤赤赤
# ２．青赤青赤赤赤
# ３．青赤赤赤赤赤

だが、以下でも確率同じ。何となく
# １．赤赤赤赤青 ★
# ２．赤赤赤赤青青
# ３．赤赤赤赤青赤
★は赤でも青でもどっちでもOK

P(# １) >P(# ２) >P(# ３) ∵青出やすい

R言語等でシミュレーションされ、
自身の確率直感が正しいのを
確認できるとは、素晴らしい。

**>>738** · 2019/12/12(木) 18:31:03.53

上記の件、若干の訂正とする

# １．赤青赤赤赤
# １'．赤赤赤赤青とすると、
# １と# １'は、直感で同じ確率と
思ってたが間違えのようだ。

当方のシミュレーションで、
# １は、0.1248
# １'は、0.1271
となった。微妙だけど、多分だ。

やっぱり確率計算をコンピュータで
モンテカルロシミュレーションのは
素晴らしい。

**１３２人目の素数さん** · 2019/12/12(木) 19:21:51.88

>>735
【grepは未使用の糞真面目な方法】

# １．赤青赤赤赤
についての確率、ほぼ厳密解を得た
# １は、0.124774…　だと思う

計算は、モンテカルロ法でない方法
でプログラム、計算した。
で、grepは使用してない。
ちなみに計算誤差は、ほぼ皆無なハズ

ソースコードイメージ
p01 = (1/3)^4*(2/3)
p02 = p01
p03 = p01 * (1 - p01)
p04 = p01 * (1 - p01 - p02)
p05 = p01 * (1 - p01 - p02 - p03)
…
p16 = p01 * (1 - p01 - p02 - p03 … - p14)
とし、
p01~p16 の合計を算出したところ、
0.124774…　　となった

**１３２人目の素数さん** · 2019/12/13(金) 07:46:08.65

>>738
レスありがとうございます。

私の直感

# １．赤青赤赤赤
# ２．青赤青赤赤赤
# ３．青赤赤赤赤赤

# １．を６個に書き換えて　#２．と並べると
# １．★ 赤青赤赤赤
# ２．青赤青赤赤赤
★は赤でも青でもどっちでもOKだから#１．の方が起こりやすい

# ２と# ３を比べると
# ２．青赤青赤赤赤
# ３．青赤赤赤赤赤
３個めでは青の方がでやすいので
ら#２．の方が起こりやすい

よって、P(# １) >P(# ２) >P(# ３)

**１３２人目の素数さん** · 2019/12/14(土) 03:40:27.74

分からない問題スレから、
>>
1回3.6%で激レアが出るガチャを10回回した確率って
36%なのでしょうか？
それとも0.964*0.964*0.964(略 0.964を10回電卓にかけた数なのでしょうか？
教えてください。
<<

百万回シミュレーション

p=3.6/100
N=10
sim <- function() any(rbinom(N,1,p)==1)
mean(replicate(1e6,sim()))

> mean(replicate(1e6,sim()))
[1] 0.306628

**１３２人目の素数さん** · 2019/12/14(土) 03:51:08.05

シミュレーション　その２

gacha=c(rep(1,36),rep(0,1000-36))
sim2 <-function() any(sample(gacha,10,replace=TRUE)==1)
mean(replicate(1e6,sim2()))

> mean(replicate(1e6,sim2()))
[1] 0.306904

何故か0.36にならない、どうしてだろ？

**１３２人目の素数さん** · 2019/12/14(土) 04:10:22.87

シミュレーションその３　（処理速度の関係で１０万回の平均）

gacha=c(rep(1,36),rep(0,1000-36))
sim3 <- function() any(replicate(10,sample(gacha,1))==1)
mean(replicate(1e5,sim3()))

> mean(replicate(1e5,sim3()))
[1] 0.30691

これも0.307弱だな。

何が悪いんだろ？　俺の頭かな？？

**１３２人目の素数さん** · 2019/12/14(土) 05:45:07.83

3.6%のガチャを10回回して全部外れる確率は
(1 - 0.036) ^ 10 ≒ 0.6930592

したがって、
3.6%のガチャを10回まわして1回以上当たる確率は
1 - 0.6930592 = 0.3069408

**１３２人目の素数さん** · 2019/12/14(土) 08:18:13.50

>>745
ありがとう、
やっぱ、
これだったｗ
何が悪いんだろ？　俺の頭かな？？

**１３２人目の素数さん** · 2019/12/14(土) 08:19:41.64

>>746
シミュレーションにバグはなかったｗ

**１３２人目の素数さん** · 2019/12/14(土) 12:29:26.79

小学生でも解けるぞ。がんばれ。

**１３２人目の素数さん** · 2019/12/14(土) 16:29:58.85

確率統計なんぞ無意味つまりインチキ
事実、結果こそすべて

**１３２人目の素数さん** · 2019/12/14(土) 19:21:15.54

その確率計算の激ナイーブな解法を示す

激レアは、レアだから1個とみなす。
故に、母集団の個数は
1÷0.036 = 27.777…　きっと28個だ

手順1) 28枚のカードがある
手順2) 重複しない1~28の番号を振る
手順3) 28枚のカードをシャッフル
手順4) 1~10枚目のどれか1となる確率
と絶対同じハズ、だから、

P(1枚目で当) =27P27 ÷ 28P28 = 1/28
P(2枚目で当) も同様に1/28
P(3枚目で当) も同様に1/28
…
P(10枚目で当) も同様に1/28

で此等10個の事象は背反事象だから、
P = 10/28 = 0.357
　∵有効数字3桁と勝手にしちゃう

ところてガチャってゲーム何か
よくわかんないけど、計算しちゃった

**１３２人目の素数さん** · 2019/12/14(土) 20:44:12.02

追記というか突然ですが、
そのガチャ3.6%の件、
超幾何分布なのか。二項分布なのか。
確率の小さいとか、母集団が小さい　とかだと無視できないと思われる。

確率統計はギャンブル派生数学ぢゃ。
生半可な知識ではカモにされる。

現代の若者たちは、数学特に
確率統計はじめとするギャンブル
の能力が特段に欠けており、
R言語等のプログラミング教育で
ギャンブルゲームを学習すべきだ。

奇麗事の学問だけの今日の数学ぢゃ
カモにされるだけ。
健全な娯楽として賭博系確率統計学
をC R Java Pyson Javascript BASIC
の何れかを学校で学習すべきだ。

**１３２人目の素数さん** · 2019/12/16(月) 09:18:31.40

私のガチャのイメージは、
かつて、昭和の駄菓子屋によくある
ガチャガチャすなわちカプセルトイ

そのような健全的なギャンブルが
超幾何分布とか二項分布の理解に
役立つのだ。
限られたお小遣の10円玉数枚で、
如何にレアアイテムをGetするかを
子供らは、思考するからだ。
「残り物には福がある」は
確率統計的には正しいのか
子供同士で文学的に議論したものだ。

さて、今のガチャは恐らくは、
デジタルの媒体のスマホゲームだ。
二項分布でよいだろう。
R言語には、
二項分布の密度関数、それの累積関数
はモチロン、それに従う乱数生成を
提供してるようだ。

ゲームの仕組みが複雑化する今、
R言語等の乱数生成プログラミングで
これからのデジタル化ギャンブル社会
で、お金より大切な激レアをドンドン
無限にゲットできる人材の増加を
期待できる可能性を秘めている

**１３２人目の素数さん** · 2019/12/17(火) 06:37:18.85

カプセルトイ（ガチャ）１台に1000個カプセルが入っていて36個がアタリ(レアアイテム）とする。
同じカプセルトイが10台ある。カプセル取り出し後は補充されない。
アタリを１個でも手に入れる確率は
１台から10個取り出す場合(G10)と
１台から１個を１０台で取り出す場合(G01)
ではどちらが高いか？

そのシミュレーション

rm(list=ls())
N=1000 ; K=36 ;n=10 # アタリ3.6％
g=rep(c(1,0),c(K,N-K))
G10 <- function() any(sample(g,n,replace=FALSE)>0)　# 非復元（超幾何分布）
G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
mean(replicate(1e6,G10()))
mean(replicate(1e6,G01()))

**１３２人目の素数さん** · 2019/12/17(火) 08:10:38.97

シミュレーションと理論値

> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.307745
[1] 0.3081121
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.307295
[1] 0.3069408

１台から１０個取り出した方がいいみたい。

**１３２人目の素数さん** · 2019/12/17(火) 14:14:51.68

>>754 なるほど、
どの台も全部で、1000カプセルで、
どの台も当りが、36カプセルだと
G10(1台で10個)取り出す方が
僅かですが、確率良さそうですね。

シミュレーションも理論値も
同様な結論のようであり、
G10(1台で10個)取り出す方が僅かに
有利が分かり何か楽しかったです。

仮に、もしどの台も
250カプセル中当たり9カプセルなら、さらにG10戦略がG01戦略より有利な
感触を掴めました。

**１３２人目の素数さん** · 2019/12/17(火) 15:44:13.73

>>755
1000個だとシミュレーションの差が微妙で再現性が不安だったけど
250個にすると、差が明らかにつきますね。

> N=250 ; K=9 ;n=10 # アタリ3.6％
> g=rep(c(1,0),c(K,N-K))
>
> G10 <- function() any(sample(g,n,replace=FALSE)>0)　# 非復元（超幾何分布）
> G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
>
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.31082
[1] 0.3117069
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.306782
[1] 0.3069408

**１３２人目の素数さん** · 2019/12/18(水) 14:18:29.16

rhyper　で超幾何分布の乱数を発生させることができたんだな。
これと　rbinomを使うと　sample関数を使わなくても同じことができた。

N=250 ; K=9 ;n=10 # アタリ3.6％

mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n) # 非復元（超幾何分布）
mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n # 復元(二項分布)

結果
> mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.311581
[1] 0.3117069
> mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n
[1] 0.307213
[1] 0.3069408

1台から順次１０個取り出すのをイメージすれば残り物には福があるとも言えなくもないな。

>

**１３２人目の素数さん** · 2019/12/18(水) 21:18:37.22

"フランスの数学者パスカル（1623～1662)が1654年にフェルマーにあてた手紙が、現在の確率
論の始まりだと言われている。当時の有名な賭博師メレがパスカルに以下のような問題を持ち
込み、その問題についてがその手紙のやりとりの中で論じられているそうである。
甲乙二人がおのおの32ピストル（当時のお金の単位）の金を賭けて勝負したとする。
そしてどちらかが先に３点を得たものを勝ちとし、勝った方がかけ金の総額64ピストルをもら
えるとする。ところが甲が２点、乙が１点を得たとき、勝負が中止になってしまった。
このとき、二人のかけ金の総額64ピストルを甲と乙にどのように分配すればよいだろうか。
ただし二人の力は互角で、勝つ確率はそれぞれ1/2ずつだとする。"

# 先にw(=3)点を得たものを勝ち、甲がA(=2)点、乙がB(=1)点を得たとき、勝負が中止
# 甲の確率を求める

gambling <- function(A=2,B=1,w=3,k=1e5){ # k : how many simulate
sim <- function(){
while(A < w & B < w){
g = rbinom(1,1,p=0.5)
if(g==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}

> gambling(2,1,3)*64 #64ピストル配分
[1] 48.13056
> gambling(2,0,4) # 日本シリーズ
[1] 0.80945

**１３２人目の素数さん** · 2019/12/19(木) 09:59:55.45

>>758の計算結果は、興味深い結果だ
で、
自己流の直感でラフな計算をしてみた。
甲1点、乙1点 ⇒1:1按分 ∵自明　★
甲3点、乙1点 ⇒1:0按分 ∵自明　☆

安直に★と☆の平均を考えると
甲2点、乙1点 ⇒1:0.5按分 ∵安直

で、総額で64ピストルだから、
甲に、64*1/1.5 = 42.67 ピストル
乙に、64*0.5/1.5 = 21.33 ピストル
を配分すれば、思ってしまった。

ラフな計算だと、少々不公平なことが
起こるようだ。

R言語等のシミュレーションは、
超天才のワィの直感より、
さらに厳密な確率計算ができそうだ。

**１３２人目の素数さん** · 2019/12/19(木) 16:21:08.12

>>758
最後はAが得点して勝負が決定するのを忘れずにAが勝者になるまでのBの点数で場合分け

NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B：現在の点数 ,p:甲の勝率
ans=0
for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数
ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A)
}
return(ans)
}

> NS(3,2,1,0.5)*64 #64ピストルの分配
[1] 48
> NS(4,1,0,0.5) #日本シリーズで先勝したほうが優勝する確率
[1] 0.65625

NSは日本シリーズの頭文字から命名、No Skinではないｗｗ

**１３２人目の素数さん** · 2019/12/19(木) 18:00:34.70

賭博等の確率計算をイメージすると
ワィの脳内は活性化される　だから、
ドンドン、妄想全開となるぅぅぅ。

という訳で、妄想内容を記載する

17世紀の地球に行き、胴元になって
賭博ビジネスでガッポリ儲けるゼェ。
17世紀の地球で、胴元になって、
コイントス賭博事業をするゼェ。

で、システムは、2回コイントスし、
表の回数が0回⇒参加者が32万円GET
表の回数が1回⇒胴元に、32万円PAY
表の回数が2回⇒参加者が32万円GET
参加料は、1万円

まんまと、参加者が沢山 1000人いた
ある参加者は、以下の様に怪説をした

組み合わせは、3通りだ。
1通り目　表の回数が0回
2通り目　表の回数が1回
3通り目　表の回数が2回
だから、参加者の勝つ確率は2/3だ

で、モチロン、胴元ワィの商売は成功
およそ、1000万円儲かっちゃた。

**１３２人目の素数さん** · 2019/12/20(金) 02:07:38.50

0745
しろ@huwa_cororon 11月27日
苦節6ヶ月、初満点&一等賞です！
https://twitter.com/huwa_cororon/status/1199593474128896000
https://twitter.com/5chan_nel (5ch newer account)

**１３２人目の素数さん** · 2019/12/21(土) 10:49:06.51

日本シリーズで賭けをする。

日本シリーズは先に4勝したチームが優勝。
勝率はそれまでの（引き分けを除いた）通算勝率に従うとする。
シリーズ開始前の通算成績はA:2勝、B:4勝であった。
今シリーズでAが先勝（第一試合に勝利）した。
この時点でA,Bのどちらに賭ける方が有利か？

**１３２人目の素数さん** · 2019/12/21(土) 22:03:14.56

分からない問題スレで３進法の小数とかが話題になっていた

https://rio2016.5ch.net/test/read.cgi/math/1567866548/845
>>
845 名前：１３２人目の素数さん[] 投稿日：2019/12/21(土) 05:52:50.74 ID:wum1jR1j
...
10進数で表された5/9を3進数になおせという問題があって、答は0.12だったのですが、
...
<<

こんなのを作ってみた。

# 小数点付きの数numをN進法で表示する

rm(list=ls())
dec2basen <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納

x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
if(N<=10){ # Nが10以下は数値として表示
cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='')
}
else{ # Nが1１以上は整数部分と小数部分を数列で表示
print(list(int=r,decimal=a))
}
invisible()

}

> dec2basen(5/9,3)
0.120000000000000

> dec2basen(3.14159265359,7)
3.06636514320
> dec2basen(3.14159265359,8)
3.11037552421
> dec2basen(3.14159265359,9)
3.12418812407
> dec2basen(3.14159265359,14)
$int
[1] 3

$decimal
[1] 1 13 10 7 5 12 13 10 8 1 4

**１３２人目の素数さん** · 2019/12/22(日) 08:18:39.30

>>764
バグの原因を追求したら、Rの仕様のせいみたいだな。

> (1.2-1)*5==1
[1] FALSE
> (1.2-1)*5>1
[1] FALSE
> (1.2-1)*5<1
[1] TRUE

これがfloor関数が誤値を返してくる原因だった。
> floor(1)
[1] 1
> floor((1.2-1)*5)
[1] 0

**１３２人目の素数さん** · 2019/12/22(日) 08:52:46.43

>>765
数値演算についてもっと勉強しな。
たいていの言語でそうなる。

**１３２人目の素数さん** · 2019/12/22(日) 13:06:36.02

すると、pythonに移植しても同じ結果ということか。

IPython 6.5.0 -- An enhanced Interactive Python.
C:\bin\Anaconda3\lib\site-packages\ipykernel\parentpoller.py:116: UserWarning: Parent poll failed. If the frontend dies,
the kernel may be left running. Please let us know
about your system (bitness, Python, etc.) at
ipython-dev@scipy.org
ipython-dev@scipy.org""")

(1.2-1)*5==1
Out[1]: False

(1.2-1)*5>1

Out[2]: False

(1.2-1)*5<1

Out[3]: True

(1.2-1)*5
Out[4]: 0.9999999999999998

**１３２人目の素数さん** · 2019/12/22(日) 17:15:43.69

round を使って回避することにした
> 0.728*5-3.64
[1] -4.440892e-16
> round(0.728*5-3.64,5)
[1] 0

んで、デバッグ版

# 小数点付きの数numをN進法で表示する(62進法まで対応 0-9,a-z,A-Z)
rm(list=ls())
dec62 <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示数列を格納
k=max(nchar(num)-nchar(floor(num))-1,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
x=round(num-floor(num),k) # e.g. 7.28-floor(7.28)-0.28 != 0に対応
for(i in 1:k){
y=round(x*N,k) # e.g. 0.728*5-3.64 !=0 に対応
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters,LETTERS)[1:N]
if(N<=62){ # Nが62以下は数値として表示
cat(paste(fig[b$integer+1],sep='',collapse=''),
'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが63以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b) # b$num:検証用
}

**１３２人目の素数さん** · 2019/12/22(日) 17:25:01.53

>>767
Haskellでもやってみた

Prelude> (1.2-1)*5
0.9999999999999998

Prelude> 0.72*5-3.6
-4.440892098500626e-16

**１３２人目の素数さん** · 2019/12/25(水) 07:12:30.44

>>766
レスありがとうございました。
アドバイスに従ってちょっと勉強してみました。
他の言語に移植しても無駄とわかって助かりました。
こうなる理由が理解できました。

(1+1/10-1)*10==1
[1] FALSE
> (1-1+1/10)*10==1
[1] TRUE

**１３２人目の素数さん** · 2020/01/06(月) 12:56:45.01

からnまでの順列を列挙するスクリプト。
これをwhile使って高速化するにはどうすればいいだろう？
俺の頭では思いつかない。

perm <- function (n) {
v=1:n
sub <- function(n, v) {
if (n == 1)
matrix(v, 1, 1)
else {
x = NULL
for (i in 1:n) x =rbind(x, cbind(v[i], sub(n -1,v[-i])))
x
}
}
sub(n, v)
}

> perm(4)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 1 2 4 3
[3,] 1 3 2 4
[4,] 1 3 4 2
[5,] 1 4 2 3
[6,] 1 4 3 2
[7,] 2 1 3 4
[8,] 2 1 4 3
[9,] 2 3 1 4
[10,] 2 3 4 1
[11,] 2 4 1 3
[12,] 2 4 3 1
[13,] 3 1 2 4
[14,] 3 1 4 2
[15,] 3 2 1 4
[16,] 3 2 4 1
[17,] 3 4 1 2
[18,] 3 4 2 1
[19,] 4 1 2 3
[20,] 4 1 3 2
[21,] 4 2 1 3
[22,] 4 2 3 1
[23,] 4 3 1 2
[24,] 4 3 2 1
>

**１３２人目の素数さん** · 2020/01/06(月) 21:34:47.52

>>771
自己解決　e1071というパッケージに再帰呼び出しなしのpermutationがありました。

library(e1071)
permu <- function (n) {
z <- matrix(1)
for (i in 2:n) {
x <- cbind(z, i)
a <- c(1:i, 1:(i - 1))
z <- matrix(0, ncol = ncol(x), nrow = i * nrow(x))
z[1:nrow(x), ] <- x
for (j in 2:i - 1) {
z[j * nrow(x) + 1:nrow(x), ] <- x[, a[1:i + j]]
}
}
return(z)
}

**１３２人目の素数さん** · 2020/01/09(木) 14:28:06.78

可変引数の扱い
Cで
#include <stdio.h>
int main(int argc, char* argv[])
{
while(argc--)
printf("%s\n", *argv++);
return 0;
}
をRでやるには ... を使って、c(...)もしくはlist(...)で引き出せばいいんだな。
検索してもなかなかわからなかった。

main <- function(...){
argv=c(...)
for(i in 1:length(argv)) cat(argv[i],'\n')
}

**１３２人目の素数さん** · 2020/01/15(水) 18:08:54.32

>>546
library(Rmpfr)

mpfr(factorial(50),5000)

one = mpfr(1, 5000)
factorial(one*50)

> mpfr(factorial(50),5000)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713018969967457666435945132957882063457991132016803840
> one = mpfr(1, 5000)
> factorial(one*50)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713378043612608166064768844377641568960512000000000000

何故か、oneを高精度に設定しておくと、正しい計算をしてくれた。

**１３２人目の素数さん** · 2020/01/26(日) 16:33:21.98

RのroundのIEEE仕様

> round(2.5)
[1] 2
> round(1.5)
[1] 2

は丸め誤差が減るという説明だったので実験してみた。

# 四捨五入 vs ＩＥＥＥ with FUN(=mean) for x(=c(-0.9,-0.8,...,0.8,0.9))
comp <- function(n=10,FUN=mean,x=sample((-9:9)/10,n,rep=T),print=T){
X=FUN(x)
dif = abs(X-FUN(round(x))) - abs(X-FUN(f45(x))) # round後の実行と四捨五入後の実行の差
if(dif!=0 & print) cat(dif<0,' : ',sort(x),'\n') # 差があれば表示 round後が小さければTRUE
return(dif) # 差を返す dif<0:round優位 dif>0:四捨五入優位
}
comp()

k=1e5
# mean
re=replicate(k,comp(print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))

# squared sum
comp(FUN=function(x) sum(x^2))
re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))

> # 平均 mean
> re=replicate(k,comp(print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ四捨五入引き分け
0.25235 0.16199 0.58566

> # 平方和 squared sum
> re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ四捨五入引き分け
0.40519 0.01413 0.58068

平均だと小差だが、統計計算の内部処理で使う平方和だと大差がついたので、体感できた。

**１３２人目の素数さん** · 2020/02/08(土) 17:42:30.00

"隔離中のクルーズ船では船内の換気が共通らしいから運悪く13日後に発症した奴がいるとその近くの部屋のやつは
プラスで14日隔離しないといけない。それが今の船内の状況という。

両隣のどちらかが感染したら14日延長、どの部屋も1日で感染する確率pは1%
部屋の配置は長方形でどの部屋にも両隣があるとする。
発症するか、隔離期間が終われば下船できる。
全員定員１の個室として客と乗務員を合わせた人数nは3000人。
クルーズ船から全員下船できる日数の期待値は？"

**１３２人目の素数さん** · 2020/02/08(土) 17:42:46.53

rm(list=ls())

p=0.001 # 1日の感染確率
n=3000 # 収容者人数
d14=14 # 隔離日数
Q0=rep(d14,n) # 必要隔離日数の配列の初期値
#
ODA <- function(X){ # One Day After 必要隔離日数の配列 x から１日後の配列を返す
Q=X # Xを作業用Qに代入
S=which(Q!=0) 　# 未感染者susceptibleのindex
s=length(S) # 未感染者人数
m=rbinom(1,s,p) #　その日の感染数
if(m==0){ # 感染者0なら未感染者の
Q[S]=Q[S]-1　 # 必要隔離日数を１日減らして返す
invisible(Q)
}else{
I.idx=sample(S,m) # 未感染者Sからm人が感染、感染者のindex (Infected.idx)
id2ext <- function(id){ # 感染者idからの両隣のidを返す
plus =ifelse(id==n,1,id+1) # 最大番号なら1を返す,それ以外は1増やす
minus=ifelse(id==1,n,id-1) # １なら最大番号を返す,それ以外は１減らす
unique(c(plus,minus)) # 重複を除く
}
E.idx=as.vector(sapply(I.idx,id2ext)) # 感染者の両隣のidを列挙し
E.idx=unique(E.idx) # 重複を除く Extended quarantine.idx
Q[E.idx]=d14 # 両隣は必要隔離日数をd14にリセット
Q[I.idx]=0 # 感染者は必要隔離日数=0 (リセット者に重複していたら０で上書き)
IE.idx=unique(c(I.idx,E.idx)) # 感染もしくは感染者の隣のindex
minus1ifnotzero <- function(x) ifelse(x==0,0,x-1) # ゼロでなければ１を減じて返す
Q[-IE.idx]=minus1ifnotzero(Q[-IE.idx]) # 感染もなく感染者の隣りでもない人の隔離日数を１日減らす
invisible(Q)
}
}
sim <- function(){
AZ=FALSE # All Zero 隔離日数がすべて０か？
i=0 # カウンター
Q=Q0 # 隔離日数配列初期値
while(!AZ){ #　隔離日数がすべて0になるまで
Q=ODA(Q)　# One Day After関数を実行
AZ=all(Q==0) # All zero？
i=i+1 # カウンター
}
return(i) # 何日かかったかを返す
}

k=1e4 # k回のシミュレーションの
RE=replicate(k,sim())
mean(RE) # 平均値＝期待値

**１３２人目の素数さん** · 2020/02/08(土) 18:25:56.61

感染確率0.1％やん
低すぎぃ

**１３２人目の素数さん** · 2020/02/08(土) 19:50:08.59

>>778
一日の感染確率。
3000人いて1日3人くらい報告されたから0.1%/日くらいかな

**１３２人目の素数さん** · 2020/02/08(土) 19:51:16.44

実行したら全員が下船できるまで約1ヶ月になった。
コードにまだバグがあるかも。

**１３２人目の素数さん** · 2020/02/09(日) 21:31:53.47

Rのcsvの文字化け治らねえ
csvそのものをutf-8にしたら今度は不正なマルチバイトがあるってエラー
どうすればいいか教えて

**１３２人目の素数さん** · 2020/02/09(日) 22:07:23.07

>>781
どんなデータ？
ヘッダーだけじゃなくてセル内も日本語？

**１３２人目の素数さん** · 2020/02/09(日) 23:10:12.62

>>781
エスパーじゃないから、まず確認だけど

・Rを使ってる環境(OS)とRのバージョンは？
・CSVのコード変換に使ってるソフトは？
・CSVを読み込んでる関数と指定している引数は？

**１３２人目の素数さん** · 2020/03/06(金) 19:53:26.32

nCov2019 for studying COVID-19 coronavirus outbreak https://guangchuangyu.github.io/nCov2019/
後で遊んでみよう

**１３２人目の素数さん** · 2020/03/07(土) 09:46:33.43

>>784
エラーがでた。R 3.6.3にアップデートしたけど
> library('remotes')
> remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
Downloading GitHub repo GuangchuangYu/nCov2019@master
Skipping 1 packages not available: BiocStyle
Installing 5 packages: downloader, cowplot, maps, magick, BiocStyle
Error: Failed to install 'nCov2019' from GitHub:
(converted from warning) package ‘BiocStyle’ is not available (for R version 3.6.3)

**１３２人目の素数さん** · 2020/03/07(土) 11:22:48.00

>>785
エラーと通り、BiocStyleパッケージがないので別途、bioconductorからインストールすればいい。

https://bioconductor.org/packages/release/bioc/html/BiocStyle.html

**１３２人目の素数さん** · 2020/03/07(土) 11:34:31.80

>>786
ありがとうございます。

**１３２人目の素数さん** · 2020/03/07(土) 11:38:15.72

We also developed a web app (http://www.bcloud.org/e/) with interactive plots and simple time-series forecasts.

という記述があるんだけど　関数 nls による非線形回帰分析　で作成しているのかなぁ？

**１３２人目の素数さん** · 2020/03/09(月) 20:05:18.91

面白い問題スレをRでシラミ潰しにカウントしてみた。

"縦n個、横n個のマス目のそれぞれに 1,2,3,...,n の数字を入れていく。
このマス目の横の並びを行といい、縦の並びを列という。
どの行にも、どの列にも、2つの対角線上にも同じ数字が1回しか現れない入れ方は何通りあるか求めよ。(2020京大文系　改)
n=4がオリジナル
"
library(gtools)
n=4
(a=permutations(n,n)) # ｎの順列
r=nrow(a)
f<-function(x){ # x=c(1,2,3,4) -> rbind(a[1],a[2],a[3],a[4])
ans=NULL
for(i in x){
ans=rbind(ans,a[i,])
}
return(ans)
}

b=permutations(r,n)
head(b) ; tail(b)
B=apply(b,1,f)
diag2 <- function(x){ # 行列xの右上から左下への対角線の配列を返す
n=nrow(x)
ans=numeric(n)
for(i in 1:n) ans[i]=x[i,n+1-i]
return(ans)
}
g <- function(v){ # matri(v,n)で列、対角線の要素がすべて異なればTRUEを返す
n=sqrt(length(v))
(x=matrix(v,n))
if(length(unique(diag (x))) < n) return(FALSE)
if(length(unique(diag2(x))) < n) return(FALSE)
flag=TRUE
for(i in 1:n){
if(length(unique(x[,i])) < n){
flag=FALSE
break
}
}
return(flag)
}
counter=NULL
for(i in 1:ncol(B)){
if(g(B[,i])) counter=append(counter,i)
}
length(counter)
counter
matrix(B[,counter[1]],n)
matrix(B[,counter[2]],n)
matrix(B[,counter[48]],n)

**１３２人目の素数さん** · 2020/03/10(火) 07:11:26.36

>>789
n=5 にするとメモリー不足でエラーがでた。
サンプリングでやってみる
# simulation
n=5
sim <- function(n=5){
v=as.vector(t(replicate(n,sample(n)))) # n × n行列要素をベクトル化
if(g(v,dia=TRUE)){ # diagonal latin squareなら
print(matrix(v,n)) # 行列化して表示
return(1)
}
return(0)
}
k=1e6
sum(replicate(k,sim(5)))

**１３２人目の素数さん** · 2020/03/10(火) 12:07:20.32

１億回モンテカルロしてようやく、1個みつかった。

[,1] [,2] [,3] [,4] [,5]
[1,] 2 5 4 3 1
[2,] 4 3 1 2 5
[3,] 1 2 5 4 3
[4,] 5 4 3 1 2
[5,] 3 1 2 5 4

**１３２人目の素数さん** · 2020/03/11(水) 11:06:12.00

嫌ね最近は総当たりアルゴリズムで枝刈り使えば実用的な時間で解決できるこ知らない人が多くて

**１３２人目の素数さん** · 2020/03/15(日) 08:45:11.95

versionが3.6,xになって、str2lang とか　str2expressionという関数が追加されたんだな。

https://id.fnshr.info/2019/07/20/r3-6-0/
で知った

例

rep("1:6",5)
paste(rep("1:6",5),collapse=',')
str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))
eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')))

(s=rep("1:6",5))
(str=paste(s,collapse=','))
(lang=str2lang(paste("expand.grid(",str,')')))
eval(lang)

**１３２人目の素数さん** · 2020/03/21(土) 17:23:30.68

The Clinical and Chest CT Features Associated with Severe and Critical COVID-19 Pneumonia

https://journals.lww.com/investigativeradiology/Abstract/publishahead/The_Clinical_and_Chest_CT_Features_Associated_with.98832.aspx

の付表にこんなデータがあった

# Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value
# Number of lobes Median (interquartile range, IQR)
# 5 (4-5) 5 (5-5) 5 (3-5) 0.003

重症群２５例　中等症群５８例でCTで肺炎像が確認できた肺葉数(１～５）が
全体で中央値が５　第一四分位数(少ない方から数えて25%の値）が４、第三四分位数(少ない方から数えて75%の値）が５
重症群では中央値５、第一四分位数５、第三四分位数５
中等症群では中央値５、第一四分位数３、第三四分位数５で
p値が０．００３であるという。
０．００２５～０．００３５までとして、

この条件に見合う８３例をリストアップすると
> re[100,][1:25]
[1] 1 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
> re[100,][26:83]
[1] 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[48] 5 5 5 5 5 5 5 5 5 5 5
という類の組み合わせで
19504通りと計算された。

重症群の肺葉数の平均は４．７６
中等症群の平均は４．０４
と計算された。

タイ値がこれだけ多いデータにWilcox検定を使う意味がわからんな。

**１３２人目の素数さん** · 2020/03/22(日) 08:57:06.28

# 感度0.7　特異度0.9のPCR検査で1000人中10人が陽性だったときの有病率の信頼区間をシミュレーションして計算するスクリプト

# アルゴリズムには自信がないｗ

sim <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){
# n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画
prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定

PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算
prev.post = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算
# NPVでも同じ結果
# NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算
　# prev.post = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算
mean=mean(prev.post) # 期待値
median=median(prev.post) # 中央値
mode=density(prev.post)$x[which.max(density(prev.post)$y)] # 最頻値
LU=unlist(HDInterval::hdi(prev.post))[1:2] # 95%信頼区間
re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値
if(print){ # ヒストグラム描画
par(mfrow=c(2,2))
hist(prev,freq=F,xlim=c(0,1)) ; lines(density(prev),lwd=2)
hist(PPV,freq=F,xlim=c(0,1)) ; lines(density(PPV),lwd=2)
hist(prev.post,freq=F,xlim=c(0,1)) ; lines(density(prev.post),lwd=2)
BEST::plotPost(prev.post)
par(mfrow=c(1,1))
}
print(round(re,4)) # 概算値表示
invisible(list(re,prev.post)) # 代表値と事後有病率を非表示で返す
}
sim(100,1)
sim(1000,10)
sim(10000,100)
sim(100000,1000)

**１３２人目の素数さん** · 2020/03/22(日) 14:13:54.45

>>795(自己レス)
間違いに気づいたので撤回します。

**１３２人目の素数さん** · 2020/03/23(月) 14:55:08.54

MCMCで解決してみた。

【富山県最強伝説】新型コロナウイルスPCR検査件数　54人　陽性0人

ある集団から54人を無作為に選んでPCR検査したら陽性0であった。
PCR検査の感度0.7　特異度0.9としてこの集団の有病率の期待値と95%信頼区間を推測せよ。

#----------------------------------------------------------------------------
# 感度SEN, 特異度SPCの検査でN人中X人が陽性であったときの推定有病率prevalence
# 弱情報事前分布はprevalence ~ dunif(0,UL) , UL:上限
#----------------------------------------------------------------------------

PCRj <- function(N,X,UL=1,SEN=0.7,SPC=0.9,verbose=TRUE){ # UL:upper limit of dunif(0,UL)
library(rjags)
modelstring=paste0('
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,',UL,')
}
')
if(verbose & UL!=1) cat(modelstring)
writeLines(modelstring,'TEMPmodel.txt')
dataList=list(n=N,x=X,sen=SEN,spc=SPC)
jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList,quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p"), n.iter=1e6, thin=10)
js=as.matrix(codaSamples)
BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) ; lines(density(js[,'prev']),col='skyblue')
round(c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]),10)
}

PCRj(100,1)
PCRj(1000,10)
PCRj(10000,100)
PCRj(100000,1000)

PCRj(54,0)
PCRj(54,0,UL=0.1)

**１３２人目の素数さん** · 2020/04/03(金) 19:48:06.63

a b c d e f g h i j
1 -0.0377710 0.02994835 0.1779426 0.29057472 0.13252796 0.3279709 -0.3534731 -0.54812874 -0.3869768 -0.1233967
2 0.2707033 0.09238235 -0.1042958 -0.08485129 -0.02590955 0.2159095 -0.3338677 0.00865952 -0.1575856 0.0953795

**１３２人目の素数さん** · 2020/04/03(金) 19:50:07.75

上のDFを

library(psych)
ICC(DF）

で級内相間求めようとすると変な値でるんだけどなんででしょう？

**１３２人目の素数さん** · 2020/04/03(金) 22:01:48.08

> DF
V1 V2
a -0.03777100 0.27070330
b 0.02994835 0.09238235
c 0.17794260 -0.10429580
d 0.29057472 -0.08485129
e 0.13252796 -0.02590955
f 0.32797090 0.21590950
g -0.35347310 -0.33386770
h -0.54812874 0.00865952
i -0.38697680 -0.15758560
j -0.12339670 0.09537950
> library(psych)
> ICC(x=DF,lmer=FALSE)
Call: ICC(x = DF, lmer = FALSE)

Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.36 2.1 9 10 0.13 -0.18 0.74
Single_random_raters ICC2 0.34 2.0 9 9 0.17 -0.26 0.74
Single_fixed_raters ICC3 0.32 2.0 9 9 0.17 -0.24 0.72
Average_raters_absolute ICC1k 0.53 2.1 9 10 0.13 -0.43 0.85
Average_random_raters ICC2k 0.51 2.0 9 9 0.17 -0.69 0.85
Average_fixed_raters ICC3k 0.49 2.0 9 9 0.17 -0.63 0.84

Number of subjects = 10 Number of Judges = 2

**１３２人目の素数さん** · 2020/04/03(金) 23:40:27.04

>>800
すいません
書き方が良くなかったです。
subject が１，２
judge a,b,c,dをきちんとfirst,second,third,fourth
と書くべきでした。
subjectとjudgeはあってると思うのです。

**１３２人目の素数さん** · 2020/04/03(金) 23:42:39.27

例えば
first second third fourth
0 0.05 0.1 0.15
0 0 0.05 0.05
0.15 0.1 0.05 0.1
0 0.05 0.2 0.05
0.1 0.15 0.1 0.05
0.15 0.05 0.1 0.05

というデータで

Call: ICC(x = over_sight)

Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.0068 1 5 18 0.43 -0.19 0.48
Single_random_raters ICC2 0.0068 1 5 15 0.44 -0.19 0.48
Single_fixed_raters ICC3 0.0068 1 5 15 0.44 -0.19 0.48
Average_raters_absolute ICC1k 0.0268 1 5 18 0.43 -1.70 0.79
Average_random_raters ICC2k 0.0268 1 5 15 0.44 -1.82 0.79
Average_fixed_raters ICC3k 0.0268 1 5 15 0.44 -1.82 0.79

Number of subjects = 6 Number of Judges = 4

**１３２人目の素数さん** · 2020/04/03(金) 23:44:30.81

となるのですが、ICC1が0.0068ってこんなに低くなるのでしょうか？
結果のとり得る値が1～0なので最大15％ぐらいの変動幅で、
制度としてはかなりいい検査方法だとおもったのですが・・・

**１３２人目の素数さん** · 2020/04/04(土) 00:07:42.72

少し理解できました。
範囲制約性の問題があるのですね。
今回のサンプルがみんな低い得点で分散が低いからICCが低く算出されていました。
ダミーデータ
1,1,1,1
を追加したらICC=0.99となりました。
得点の高いサンプルを追加しないとICCで信頼性は測れないのですね。

**１３２人目の素数さん** · 2020/04/04(土) 00:09:17.81

いま有病者と健常者の両方を計測する検査法で健常者のみのデータでICCを算出しようとしているのが間違いということですね。

上記のような場合で、健常者のみでも信頼性を算出する方法っていうのはないのでしょうか？

**１３２人目の素数さん** · 2020/04/04(土) 21:44:33.98

>>804
DF=rbind(
c(0,0.05,0.1,0.15),
c(0,0,0.05,0.05),
c(0.15,0.1,0.05,0.1),
c(0,0.05,0.2,0.05),
c(0.1,0.15,0.1,0.05),
c(0.15,0.05,0.1,0.05),
c(0.0001,0.0001,0.0001,0.0001))

とダミーに0.0001を加えても
> DF2ICC1(DF)
ICC(1,1) ICC(1,k)
0.2462707 0.5665262

> psych::ICC(DF)
type ICC F df1 df2 p lower bound
Single_raters_absolute ICC1 0.25 2.3 6 21 0.072 -0.027
Single_random_raters ICC2 0.25 2.3 6 18 0.079 -0.028
Single_fixed_raters ICC3 0.25 2.3 6 18 0.079 -0.034
Average_raters_absolute ICC1k 0.57 2.3 6 21 0.072 -0.115
Average_random_raters ICC2k 0.57 2.3 6 18 0.079 -0.122
Average_fixed_raters ICC3k 0.57 2.3 6 18 0.079 -0.154

となるから、高値データでなくてもいいみたい。

**１３２人目の素数さん** · 2020/04/24(金) 21:32:48.82

R 4.0.0 来た

**１３２人目の素数さん** · 2020/04/25(土) 06:37:54.89

SIGNIFICANT USER-VISIBLE CHANGES

https://cran.r-project.org/doc/manuals/r-devel/NEWS.html

**１３２人目の素数さん** · 2020/04/25(土) 07:20:54.46

パッケージを一括インストール

install.packages(available.packages()[,1])

**１３２人目の素数さん** · 2020/04/25(土) 07:29:33.43

>>809
これやってからでないとパッケージのコンパイルができない。

After installation is complete, you need to perform one more step to be able to compile R packages:
you need to put the location of the Rtools make utilities (bash, make, etc) on the PATH.
The easiest way to do so is create a text file .Renviron in your Documents folder which contains the following line:

PATH="${RTOOLS40_HOME}\usr\bin;${PATH}"

You can do this with a text editor, or you can even do it from R like so:

writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")

**１３２人目の素数さん** · 2020/04/25(土) 08:56:31.54

>>809
libraryのフォルダで　dir /B とやって既にinstall済のリストを作って
install.packages("rstan","ggplot2",....)とやった方がいいな。
5000ほどパッケージがあるらしいけど、全部インストールしたらどれくらいの容量が必要なんだろう？

**１３２人目の素数さん** · 2020/04/25(土) 12:31:15.73

やってみてくれｗ

**１３２人目の素数さん** · 2020/04/26(日) 21:05:35.88

>>811　これでいいかな？
targetPackages <- available.packages()[,"Package"]
newPackages <- targetPackages[!(targetPackages %in% installed.packages()[,"Package"])]
install.packages(newPackages, repos = "http://cran.r-project.org";)

**１３２人目の素数さん** · 2020/04/26(日) 22:22:15.60

高血圧で不安だという老人の経過観察入院1件。発熱してなくてほっとした。
ステーキハウスがテイクアウトを始めたからインセンティブはその足しにしようっと。

**１３２人目の素数さん** · 2020/05/03(日) 15:15:43.67

>>812
コンパイルなしでやってみた。

サイズ：19.6 GB (21,082,530,395 バイト)

ディスク上のサイズ：20.3 GB (21,805,858,816 バイト)

ファイル数: 527,982、フォルダー数: 118,521

**１３２人目の素数さん** · 2020/05/03(日) 15:37:44.00

Ubuntu 20.04でRcmdrがビルドできないんですけど、これはソフトウェアが20.04に対応するまで待つしかないですか？

**１３２人目の素数さん** · 2020/05/03(日) 20:21:16.34

>>816
自己解決しました。ライブラリの不足でした。

**１３２人目の素数さん** · 2020/05/04(月) 17:06:14.69

4.0でまたsjis関連の問題が出まくりそう

**１３２人目の素数さん** · 2020/05/07(木) 01:08:38.98

4.0になってから
col=番号指定での色が目に優しくなったみたい。

hist(runif(100),col=2)
hist(runif(100),col=rgb(1,0,0,1))
hist(runif(100),col='red')

**１３２人目の素数さん** · 2020/05/07(木) 18:12:52.58

rgbで256色でそんなに変化ないんじゃない。

**１３２人目の素数さん** · 2020/05/07(木) 19:16:11.76

>>820
colours()　と入力すると　657色から選べるから256以上ありそう。

**１３２人目の素数さん** · 2020/05/07(木) 23:34:33.14

4.0からパレット変えたんじゃなかったか？

**１３２人目の素数さん** · 2020/05/07(木) 23:39:55.93

リリースノートより

The palette() function has a new default set of colours (which are less saturated and have better accessibility properties).
There are also some new built-in palettes, which are listed by the new palette.pals() function.
These include the old default palette under the name "R3". Finally, the new palette.colors() function allows a subset of colours to be selected from any of the built-in palettes.

**１３２人目の素数さん** · 2020/05/07(木) 23:46:21.45

具体的な変更点はこのあたりを参照方
https://developer.r-project.org/Blog/public/2019/11/21/a-new-palette-for-r/index.html

**１３２人目の素数さん** · 2020/05/08(金) 00:23:14.39

binomにある信頼区間をグラフ表示させてみた。

binom.ci <- function(x,n,cl=0.95){
PEci=binom::binom.confint(x,n,conf.level=cl)
y=nrow(PEci):1
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mar=c(3,3,2,2))
par(oma=c(0,4,0,0))
plot(PEci$upper,y,xlim=c(min(PEci$lower),max(PEci$upper)),type="n",
yaxt="n",ylab="",xlab="Confidence Interval",
main=paste("C.I. for", x,"out of",n))
points(PEci$upper,y,pch="|") ; points(PEci$lower,y,pch="|")
segments(PEci$lower,y,PEci$upper,y,lwd=3,
col=sample(colours(),1))
abline(v=PEci$mean[1],lty=3)
axis(2, at=y, labels=PEci$method,as.vector,las=2,cex=0.75)
PEci
}
binom.ci(3,312)

**１３２人目の素数さん** · 2020/05/08(金) 00:48:23.14

>>824
情報、ありがとうございます。

?paletteのヘルプにあるスクリプトを実行すると体感できました。

## Demonstrate the colors 1:8 in different palettes using a custom matplot()
sinplot <- function(main=NULL) {
x <- outer(
seq(-pi, pi, length.out = 50),
seq(0, pi, length.out = 8),
function(x, y) sin(x - y)
)
matplot(x, type = "l", lwd = 4, lty = 1, col = 1:8, ylab = "", main=main)
}
sinplot("default palette")

palette("R3"); sinplot("R3")
palette("Okabe-Ito"); sinplot("Okabe-Ito")
palette("Tableau") ; sinplot("Tableau")
palette("default") # reset

**１３２人目の素数さん** · 2020/05/09(土) 20:51:55.34

rjagsのコードに

theta[i,j] <- min(1, exp(-alpha[i]*t[j]) + beta[i])といいれて、確率を１以下にして二項分布させてようとすると

Node inconsistent with parents　というエラーがでた。　stackoverflow.comも明確な答が見いだせず、あれこれいじって

どうも確率1や0で二項分布する可能性があるとエラーがでるみたい、というのを発見。

1を0.99999とするとエラーが回避できた。

1を　1 - 1e-16としてもエラー回避できた。

**１３２人目の素数さん** · 2020/05/18(月) 08:52:51.27

mathematicaは内部どうなってんだか知らないが第何種ベッセル関数とか入ってくるめちゃくちゃ重い積分を計算してくれた

**１３２人目の素数さん** · 2020/05/25(月) 12:44:17.51

"
ある人物Dが新型コロナ肺炎に罹患したとする。
行動調査によって発症前にキャバクラに行っており
接客したキャバ嬢Cが人物D発症の２日後に発症していたことがわかった。
キャバ嬢Cは人物Dから移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢Cから移された可能性もあると主張して
その確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
"
潜伏期間が従う分布はこれ
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612

Gt <- function(delay){
C=rlnorm(1e6,ln_par1,ln_par2)
D=rlnorm(1e6,ln_par1,ln_par2)
mean(C-D > delay)
}
Gt(2)

**１３２人目の素数さん** · 2020/05/25(月) 12:48:45.94

怖い本で紹介されていたパッケージ　polspline　は秀逸。
発生させた乱数から確率密度関数を作ってくれる。
複雑な分布だと収束できず作れないこともあるけど

対数正規分布の差の分布を乱数発生させてpdfが作れた。

library(polspline)
delay=2
ln_par1 = 1.434065
ln_par2 = 0.6612
c=rlnorm(1e5,ln_par1,ln_par2)
d=rlnorm(1e5,ln_par1,ln_par2)
hist(c-d, freq=F,col='white',breaks=100,ylim=c(0,0.11),
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
fit=logspline(c-d)
curve(dlogspline(x,fit),-30,30,ann=F,bty='n',add=T)
segments(delay,0,delay,dlogspline(delay,fit),pch=19,col=4)
curve(dlogspline(x,fit),delay,30,add=T,type='h',col=4)
curve(plogspline(x,fit),-30,30,lty=3,bty='n')
1-plogspline(delay,fit)
fn <- function(delay) 1- plogspline(delay,fit)
curve(fn(x),0,14, bty='n' ,xlab='Delay', ylab='Probability')

**１３２人目の素数さん** · 2020/06/01(月) 16:49:09.39

今気づいたけどdplyr 1.0.0がcranに来てた

**１３２人目の素数さん** · 2020/07/25(土) 08:42:04.47

Rのグラフの出力をGIF動画にしたかったので検索したら、このページが役立った。

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/36.html

mytest <- function() {
ANS <- readline("1+2? ")
if (ANS == "3") cat("Correct!\n")
else cat("Wrong...\n")
}
mytest()

mymessage <- function() {
Sys.sleep(3)
cat("Hello.\n")
}
mymessage()

尚、GIF作製には　
https://www.screentogif.com/
がお手軽だった。

**１３２人目の素数さん** · 2020/07/26(日) 10:34:51.02

他人のコードを読んでいたら、

'%&%' <- function(x,y) paste0(x,y)

という関数が定義されていた。

これは便利と思った。

> '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi
[1] "√2 =1.4142135623731, π = 3.14159265358979"

**１３２人目の素数さん** · 2020/07/27(月) 13:23:47.24

文字列をパイプでつないでいくのは面白いかも

**１３２人目の素数さん** · 2020/07/28(火) 04:21:24.42

'%.%' <- function(x,y) sum(x,y)
でベクトルの内積とかも

**１３２人目の素数さん** · 2020/07/28(火) 04:22:11.80

訂正
'%.%' <- function(x,y) sum(x*y)

**１３２人目の素数さん** · 2020/07/29(水) 12:22:51.45

面白いけど使いどころが難しそう
%>%パイプとの親和性が無いのが痛い

**１３２人目の素数さん** · 2020/07/29(水) 13:52:08.22

%文字%は二項演算子だから、定義次第じゃない？

**１３２人目の素数さん** · 2020/07/29(水) 18:45:33.70

tidyverseや%>%を使わないならまあ

**１３２人目の素数さん** · 2020/07/30(木) 18:02:20.97

>>833
それって意味があるの？
> paste0('√2 =', sqrt(2), ', π = ', pi)
[1] "√2 =1.4142135623731, π = 3.14159265358979"
普通にpaste0のまま使った方がタイプ量が圧倒的に少ないが。

**１３２人目の素数さん** · 2020/07/30(木) 20:10:10.19

>>840
括弧の対応で混乱しそうなときにちょっと役にたつかも。

こんな感じ

eval(str2lang("expand.grid(" %&% paste0(rep("1:2",6),collapse=',') %&% ")"))

**１３２人目の素数さん** · 2020/07/31(金) 07:51:08.50

>>841
うーん、その例だとイマイチに思える

rerun(6, 1:2) %>% expand.grid
で済むものをわざわざ複雑にしているだけに感じる

**１３２人目の素数さん** · 2020/07/31(金) 13:51:36.15

>>842
無理やり、例示した感は否めないな。
標準装備でなら、
expand.grid(replicate(6,1:2,simplify = FALSE))
ですむし。

**１３２人目の素数さん** · 2020/07/31(金) 23:36:09.03

paste()やstr_c()以外を使うなら、sprintf()かglue::glue()が候補か
特にカンマが文字列に含まれるときにはpaste()より見やすくなるはず
sprintf("√2 = %s, π = %s", sqrt(2), pi)
glue("√2 = {sqrt(2)}, π = {pi}")

**１３２人目の素数さん** · 2020/08/01(土) 18:51:48.81

>>844
確かに便利ですね。
カンマを　','で入力するのは混乱のもとでしたから。
有用な情報、ありがとうございます。'

> glue::glue("√2 = {sqrt(2)}, sin(π/3) = {sin(pi/3)}")
√2 = 1.4142135623731, sin(π/3) = 0.866025403784439

**１３２人目の素数さん** · 2020/08/06(木) 00:01:24.75

初心者質問で申し訳ありません。
ある関数を作りたいのですが、その関数に入れる数値の数はランダムな場合はどのように書けばいいのでしょうか。
例えば、入れた要素をすべて掛け合わせる関数が作りたい場合、そこに入れる要素が（２，３，４）でも、（２，３，４，５，６）でも反応するような関数の書き方を教えてくださると助かります。

**１３２人目の素数さん** · 2020/08/06(木) 00:05:38.12

ベクトルで渡せばいいのでは？
そういう話じゃない？

**１３２人目の素数さん** · 2020/08/06(木) 01:18:04.99

874さん。
847さん、回答ありがとうございます。
確かにベクトルで何かするのはわかるのですが、その部分がわからないのです。
出来れば下に問題を書くので、コードを書いていただけると非常に助かります。

問、正の整数であるベクトルを入力すると、そのベクトルの要素をすべて掛け合わせた値を結果として返す関数を作りなさい。

**１３２人目の素数さん** · 2020/08/06(木) 02:03:47.23

>>848
問題文これだけ？他に条件はないの？
標準の関数にあるけどいいのかな？

関数fを作るとして、一番いいのは標準の関数に余計な変更を加えないことだから、これが答えだけど…
f <- prod

**１３２人目の素数さん** · 2020/08/06(木) 02:11:09.57

849さん、返信ありがとうございます。

実は問２がありまして、たぶんそれは既存の関数にはないと思うので、そちらも教えていただけると非常に助かります。

問２、要素がすべて整数であるベクトルを入力すると、そのベクトルの要素の絶対値をすべて掛け合わせた値を結果として返す関数を作成せよ。

**１３２人目の素数さん** · 2020/08/06(木) 05:13:12.66

>>850
f = function (x) abs(prod(x))

**１３２人目の素数さん** · 2020/08/06(木) 05:15:08.83

>>849
while やfor もしくは再帰関数を使って作らせる問題かな？

**１３２人目の素数さん** · 2020/08/06(木) 05:27:00.71

>>846
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans
}
> f(2,3,4)
[1] 24
> f(2,3,4,5,6)
[1] 720

**１３２人目の素数さん** · 2020/08/06(木) 05:32:55.68

>>850
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans=ifelse(ans>0,ans,-ans)
ans
}

> f(-1,2,3)
[1] 6
> f(1 -2,-3,-4,-5)
[1] 60

**１３２人目の素数さん** · 2020/08/06(木) 05:34:50.08

>>850
入力はベクトルだったから、

> f <- function(v){
+ ans=1
+ for(i in v) ans=ans*i
+ ans=ifelse(ans>0,ans,-ans)
+ ans
+ }
> f(c(-1,2,3))
[1] 6
> f(c(1 -2,-3,-4,-5))
[1] 60