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

1132人目の素数さん2017/08/03(木) 19:23:12.67ID:Hq1blL0O
R は統計計算とグラフィックスのための言語・環境です。
統計計算で重宝するデータ型や、複数要素を処理する演算や関数、
解析結果を表示するグラフィックなど、多彩な機能を提供します。

●関連サイト
The R Project
http://www.r-project.org/
RjpWiki
http://www.okada.jp.org/RWiki/
リンク集
http://www.okada.jp.org/RWiki/?%A5%EA%A5%F3%A5%AF%BD%B8
※前スレ
【R言語】統計解析フリーソフトR 第5章【GNU R】
http://rio2016.2ch.net/test/read.cgi/math/1380168442/

2◆2VB8wsVUoo 2017/08/03(木) 19:24:54.74ID:qtUXVqva
★★★忖度と処世術に汚染された日本人:権威主義的な支配と損したくない人達★★★
  〜〜〜芳雄氏が言う『研究者としての基本的態度』とは一体何だろうか〜〜〜

佐藤幹夫:自分自身の素朴な疑問に真剣に耳を傾ける。⇒不滅の金字塔を打ち立てる。
糞父芳雄:人間関係を駆使し他人を操り根回しを行う。⇒ハリボテお教授として君臨。

隠蔽の財務省、嘘吐きの文科省、そして問答無用に屈服させる官邸。コレでも先進国?

(佐藤師がしてたのは本物の研究だ。だが)芳雄氏がしてたのはケケケ、ケンキュウ。
外見を繕って偉そう見せさえすれば何でもヨロシ。ほんで教授になりさえすれば研究の
中身なんて何でもヨロシ。そもそも論文なんてモンは、外国の権威ある雑誌に掲載され
さえすれば、その中身のギロンなんて何でもヨロシ。そやし適当に書いてしまえ〜〜〜
中身がダメだと知ってて、ソレでもSTAP論文を外国に投稿して受理される。発覚したら
適当に言い逃れる醜い態度。オツムのダメな大学院生に「虚偽の良品ラベル」を貼って
世間に出荷するハリボテ大学は詐欺行為そのもの。世間に媚びを売って客商売に徹し、
『売れさえすれば学生の脳の質なんて何でもヨロシ』と居直る大学。そしてブランド名
だけを見て仕入れる世間。●●は一流大学やさかい、きっと優秀なエリートやろwww

中身を何も説明しないで、問答無用に上から押し付ける。ソレをイチャモンで騒いで、
そして邪魔して潰そうとする周囲の下々。大学教員も国会議事堂も、そして馬鹿板人の
遣ってる事も皆同じだ。日本人はバカ民族であり、今は外国にもちゃんとバレてるので
海外からも軽蔑されるだけであり、そのうちにどの国からも信用されなくなるだろう。

近視眼的で打算的な人生観を息子に押し付ける父親と、大脳に栄養が足りてない連中が
跋扈する永田町や霞が関に支配される国に住む不幸、一体どうしてくれるというのか。

■■■馬鹿板をスルと安倍晋三みたいな偽善者になります。そやし止めなさい。■■■


3132人目の素数さん2017/08/04(金) 03:31:08.39ID:b3+Cy8qY
板違う

4◆2VB8wsVUoo 2017/08/04(金) 06:12:55.77ID:u5H9jwGc

5◆2VB8wsVUoo 2017/08/04(金) 16:53:10.14ID:u5H9jwGc

6◆2VB8wsVUoo 2017/08/04(金) 16:53:27.29ID:u5H9jwGc

7◆2VB8wsVUoo 2017/08/04(金) 16:53:46.10ID:u5H9jwGc

8◆2VB8wsVUoo 2017/08/04(金) 16:54:03.40ID:u5H9jwGc

9◆2VB8wsVUoo 2017/08/04(金) 16:54:19.87ID:u5H9jwGc

10◆2VB8wsVUoo 2017/08/04(金) 16:54:36.75ID:u5H9jwGc

11◆2VB8wsVUoo 2017/08/04(金) 16:54:56.00ID:u5H9jwGc

12◆2VB8wsVUoo 2017/08/04(金) 16:55:13.77ID:u5H9jwGc

13◆2VB8wsVUoo 2017/08/04(金) 16:55:35.11ID:u5H9jwGc

14132人目の素数さん2017/08/06(日) 18:24:07.34ID:oDKJI1vJ
耳栓をしたら世界が変わってワロタ

15◆2VB8wsVUoo 2017/08/06(日) 18:24:24.98ID:+CYdGQny

16132人目の素数さん2017/08/06(日) 18:25:44.67ID:oDKJI1vJ
耳栓をしたら世界が変わってワロタ

17◆2VB8wsVUoo 2017/08/06(日) 18:26:01.57ID:+CYdGQny

18132人目の素数さん2017/08/06(日) 18:29:29.99ID:oDKJI1vJ
耳栓をしたら世界が変わってワロタ

19◆2VB8wsVUoo 2017/08/06(日) 18:29:52.60ID:+CYdGQny

20132人目の素数さん2017/08/06(日) 20:43:04.35ID:oDKJI1vJ
耳栓をしたら世界が変わってワロタ

21◆2VB8wsVUoo 2017/08/06(日) 22:03:22.61ID:+CYdGQny

22132人目の素数さん2017/08/07(月) 08:49:27.85ID:0YzkEl/p
耳栓をしたら世界が変わってワロタ

23132人目の素数さん2017/08/13(日) 06:51:19.01ID:+zIOnWE6
アンドロイドにRをインストールしてみた。
テキストベースなら簡易統計電卓として使える。
依拠したのはこれ。
https://stackoverflow.com/questions/36968411/installing-r-on-android
これくらいの簡単な計算ならタブレットでできるようになった。
いちいちパソコンを立ち上げなくてもすむ。

> binom.test(420,1000)$conf.int[1:2]
[1] 0.3891836 0.4512888

24132人目の素数さん2017/08/13(日) 06:52:35.49ID:+zIOnWE6
>>297
クリップボード経由でスクリプトを実行してみたが
日本語が混ざるとエラーになるな。
フィーリングカップル5対5のカップル数とカップル成立確率の
シミュレーションを走らせてみた。

http://imagizer.imageshack.com/img924/4456/siP7dz.png

きちんと動く。

25132人目の素数さん2017/08/16(水) 21:19:26.05ID:lgZzPGKw
Unfortunately installing Debian package build-essential runs into problems (a form of a dependency hell, that is),
and installing new R packages that would require compiling C or Fortran code is not currently possible. Otherwise the installation seems to work rather nicely.
ということでパッケージのインストールはエラーになった。

26132人目の素数さん2017/08/17(木) 01:51:31.67ID:qXaIFkoQ
>>25
>dependency hell
具体性がない。どこ情報?

27132人目の素数さん2017/08/18(金) 14:16:09.61ID:1KawOE0O

28132人目の素数さん2017/08/18(金) 14:17:55.45ID:1KawOE0O
Console R はコンパイラーが使えたからPremium版を買ったけど
今は公開中止で入手できない。

29132人目の素数さん2017/08/18(金) 19:59:31.55ID:x/xXagzJ
>>27
Debianの話ではなくて、Andoridか、紛らわしい。
ってか、クロスコンパイル環境でビルドするのではなくて、
DebianのバイナリをAndoridの中のエミュレータで動かすという話?

30132人目の素数さん2017/08/18(金) 23:41:48.42ID:yku0iUDJ
こいつらスマホからR使ってんの?

31132人目の素数さん2017/08/19(土) 00:56:43.44ID:joMJEhOz
FireHDで使いたい

32◆2VB8wsVUoo 2017/08/19(土) 02:20:10.86ID:LB3Hl+jp

33◆2VB8wsVUoo 2017/08/19(土) 02:20:26.96ID:LB3Hl+jp

34◆2VB8wsVUoo 2017/08/19(土) 02:20:43.33ID:LB3Hl+jp

35◆2VB8wsVUoo 2017/08/19(土) 02:21:17.96ID:LB3Hl+jp

36◆2VB8wsVUoo 2017/08/19(土) 02:21:33.60ID:LB3Hl+jp

37◆2VB8wsVUoo 2017/08/19(土) 02:21:49.32ID:LB3Hl+jp

38◆2VB8wsVUoo 2017/08/19(土) 02:22:05.59ID:LB3Hl+jp

39◆2VB8wsVUoo 2017/08/19(土) 02:22:21.98ID:LB3Hl+jp

40◆2VB8wsVUoo 2017/08/19(土) 02:22:39.87ID:LB3Hl+jp

41◆2VB8wsVUoo 2017/08/19(土) 02:22:56.35ID:LB3Hl+jp

42132人目の素数さん2017/08/19(土) 07:33:24.79ID:xFz3lrjm
>>31
先ず、Google Playをいれて>23の手順で使えるようになるはず。

43132人目の素数さん2017/08/19(土) 07:36:42.75ID:xFz3lrjm
>>29
GNURoot Debian provides a method for you to install and use Debian and its associated applications/packages alongside Android.
とあるからエミュレータなのだろう。

44132人目の素数さん2017/08/19(土) 09:11:59.25ID:xFz3lrjm
>>31
他の泥タブレットと同じくFireHDでもできたよ。

45132人目の素数さん2017/08/21(月) 01:18:38.70ID:Dj1b7ycp
気づいたら前スレ埋められてたのか
いつも思うがRに親でも殺されたのかと

46◆2VB8wsVUoo 2017/08/21(月) 01:19:09.87ID:zqnQwMkT

47◆2VB8wsVUoo 2017/08/21(月) 01:19:28.11ID:zqnQwMkT

48◆2VB8wsVUoo 2017/08/21(月) 01:19:45.27ID:zqnQwMkT

49◆2VB8wsVUoo 2017/08/21(月) 01:20:02.14ID:zqnQwMkT

50◆2VB8wsVUoo 2017/08/21(月) 01:20:19.41ID:zqnQwMkT

51◆2VB8wsVUoo 2017/08/21(月) 01:20:38.03ID:zqnQwMkT

52◆2VB8wsVUoo 2017/08/21(月) 01:20:56.12ID:zqnQwMkT

53◆2VB8wsVUoo 2017/08/21(月) 01:21:13.27ID:zqnQwMkT

54◆2VB8wsVUoo 2017/08/21(月) 01:21:28.95ID:zqnQwMkT

55◆2VB8wsVUoo 2017/08/21(月) 01:21:46.90ID:zqnQwMkT

56132人目の素数さん2017/08/27(日) 21:42:47.08ID:RuVE0xlS
2^127-1 =170141183460469231731687303715884105727は素数

Mathematicaは正しく計算してくれるけど

Rだと

> 2^127-1
[1] 170141183460469231738248242066868668888

と表示された。

あまり大きな数字は扱えないようだ。

57132人目の素数さん2017/08/27(日) 21:46:56.21ID:RuVE0xlS
2^1279-1も素数

10407932194664399081925240327364085538615262247266704805319112350403608059673360
29801223944173232418484242161395428100779138356624832346490813990660567732076292
41295093892203457731833496615835504729594205476898112116936771475484788669625013
84438260291732348885311160828538416585028255604666224831890918801847068222203140
521026698435488732958028878050869736186900714720710555703168729087

58132人目の素数さん2017/08/27(日) 22:03:43.58ID:dgHzAQs5
なんかパッケージないの?

59132人目の素数さん2017/08/28(月) 20:44:54.60ID:waEfFbH1
gmpというパッケージがあったので使ってみた。

メルセンヌ素数算出のスクリプトを組んでみた
2^n-1 で n < 10000のとき

is.prime <- function(n){
i=2
while(i <= ceiling(sqrt(n)) && n%%i){
i=i+1
}
if(n>i && !n%%i) return(FALSE)
else return(TRUE)
}

library('gmp')

primes <- function(n){
xx=c(2,seq(3, n, by=2))
xx[sapply(xx,is.prime)]
}
N=10000
p = primes(N)
M = as.bigz(2)^p - 1
p[isprime(M)>0]

>p[isprime(M)>0]
[1] 2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279
[16] 2203 2281 3217 4253 4423 9689 9941

60132人目の素数さん2017/08/29(火) 00:26:34.49ID:hc5LENHz
Rって何のための言語だっけ

61132人目の素数さん2017/08/29(火) 00:33:01.32ID:H3vsxGlL
Sのフリー版

62132人目の素数さん2017/08/29(火) 02:25:16.29ID:Dd6KiNrk
Rに親を殺されたマンはよ

63◆2VB8wsVUoo 2017/08/29(火) 04:56:27.80ID:TbkIY/Vo

64◆2VB8wsVUoo 2017/08/29(火) 04:56:46.76ID:TbkIY/Vo

65◆2VB8wsVUoo 2017/08/29(火) 04:57:05.09ID:TbkIY/Vo

66◆2VB8wsVUoo 2017/08/29(火) 04:57:24.20ID:TbkIY/Vo

67◆2VB8wsVUoo 2017/08/29(火) 04:57:43.10ID:TbkIY/Vo

68◆2VB8wsVUoo 2017/08/29(火) 04:58:02.49ID:TbkIY/Vo

69◆2VB8wsVUoo 2017/08/29(火) 04:58:23.65ID:TbkIY/Vo

70◆2VB8wsVUoo 2017/08/29(火) 04:58:44.27ID:TbkIY/Vo

71◆2VB8wsVUoo 2017/08/29(火) 04:59:03.74ID:TbkIY/Vo

72◆2VB8wsVUoo 2017/08/29(火) 04:59:22.59ID:TbkIY/Vo

73132人目の素数さん2017/08/29(火) 11:36:27.20ID:uQrlExl4
>>60
貧乏人のためのデータプログラミング言語

74◆2VB8wsVUoo 2017/08/29(火) 12:03:33.07ID:TbkIY/Vo

75◆2VB8wsVUoo 2017/08/30(水) 04:33:42.82ID:xZ8twSlP

76◆2VB8wsVUoo 2017/08/30(水) 04:34:00.30ID:xZ8twSlP

77◆2VB8wsVUoo 2017/08/30(水) 04:34:17.96ID:xZ8twSlP

78◆2VB8wsVUoo 2017/08/30(水) 04:34:35.92ID:xZ8twSlP

79◆2VB8wsVUoo 2017/08/30(水) 04:34:53.84ID:xZ8twSlP

80◆2VB8wsVUoo 2017/08/30(水) 04:35:11.78ID:xZ8twSlP

81◆2VB8wsVUoo 2017/08/30(水) 04:35:29.01ID:xZ8twSlP

82◆2VB8wsVUoo 2017/08/30(水) 04:35:46.61ID:xZ8twSlP

83◆2VB8wsVUoo 2017/08/30(水) 04:36:05.83ID:xZ8twSlP

84132人目の素数さん2017/08/30(水) 21:00:36.30ID:XCGWMXRG
荒らし(◆2VB8wsVUoo)の正体は元筑波大学准教授で数学者の増田哲也。
増田哲也は2007年に痴漢で逮捕され、精神を病んで2ch数学板を荒らすようになった。
自ら増田哲也とカミングアウトしている。父は植物学者の増田芳雄。
荒らしが酷く、数学板で専用スレが10スレ以上立てられた。
↓確認できる最初のスレ
http://science6.2ch.net/test/read.cgi/math/1243605006/


 ■徳島で痴漢の准教授を解雇 筑波大 (2007年8月5日 毎日新聞)
徳島県警阿南署などは5日未明、東京都足立区千住寿町、筑波大学准教授、増田哲也容疑者(50)を
県迷惑行為防止条例違反(痴漢行為)容疑で逮捕した。

調べでは、増田容疑者は、4日午後4時20分ごろから約50分にわたり、JR牟岐線の列車内で、県内の
専門学校生の女性(21)の胸や太ももなどを触った疑い。調べに対し、「夏休み期間に、講演活動を兼ね
て旅行していた。好みの女性だったのでムラムラした」と話しているという。
http://science6.2ch.net/test/read.cgi/math/1262704792/


もともとは「狢」や「猫」、「狸」などと名乗っていた。トリップも変わっている。
http://wc2014.2ch.net/test/read.cgi/math/1384592236/
増田哲也→猫◆→狢◆→狸◆


☆数学コテ紹介☆

・猫
本名、増○哲也。筑波大准教授の頃、徳島で痴漢をやらかし職を失う。
それを期に数学界から身を引いた(←アフォー)。
逮捕されてからは精神を病み、2ch荒らしが生きがいとなった。また父の芳雄に虐待されたと思い込んでいる。
もともとは大数学者アラン・コンヌに直々教えを乞うていたほどのやり手(らしい)。
数々のお涙頂戴昔話には定評あり。本人曰く「しつこさ」のみが自身の売りなのだそう。
馬鹿を煽って2ch潰しをしているのだそう。規制されてもプロバイダーを変えて復活する。
最近、院生disが激しい。日本語、英語、フランス語をしゃべる。だが、猫語はしゃべらない。
率直に言って、客観的に人生を無駄に過ごしている希ガスる。
2ch潰しどころか数学板さえ潰すことの出来ないでいる希ガスる。
作戦倒れしている希ガスる。
でも本人はなんとも思ってないんだろうな。哀れ。哀れ。哀れ。

85132人目の素数さん2017/09/02(土) 12:47:15.63ID:3V8qFOPU
耳栓をしたら世界が変わってワロタ

86◆2VB8wsVUoo 2017/09/02(土) 13:20:51.69ID:z17/uuYO

87◆2VB8wsVUoo 2017/09/02(土) 15:52:21.39ID:z17/uuYO

88◆2VB8wsVUoo 2017/09/02(土) 19:54:32.55ID:MGAE+mBM

89◆2VB8wsVUoo 2017/09/02(土) 23:31:02.19ID:MGAE+mBM

90◆2VB8wsVUoo 2017/09/03(日) 11:06:09.04ID:TZkGcAET

91◆2VB8wsVUoo 2017/09/03(日) 14:15:50.87ID:TZkGcAET

92◆2VB8wsVUoo 2017/09/03(日) 20:37:27.93ID:TZkGcAET

93◆2VB8wsVUoo 2017/09/04(月) 02:38:55.34ID:xP4OelQr

94132人目の素数さん2017/09/04(月) 09:59:51.49ID:ZYM7kfu0
耳栓をしたら世界が変わってワロタ

95◆2VB8wsVUoo 2017/09/04(月) 10:34:11.34ID:xP4OelQr

96132人目の素数さん2017/09/04(月) 10:42:05.51ID:ZYM7kfu0
耳栓をしたら世界が変わってワロタ

97◆2VB8wsVUoo 2017/09/04(月) 10:42:22.86ID:xP4OelQr

98132人目の素数さん2017/09/04(月) 10:48:55.99ID:ZYM7kfu0
耳栓をしたら世界が変わってワロタ

99◆2VB8wsVUoo 2017/09/04(月) 10:49:35.67ID:xP4OelQr

100◆2VB8wsVUoo 2017/09/04(月) 12:08:42.71ID:xP4OelQr

101◆2VB8wsVUoo 2017/09/04(月) 12:08:59.34ID:xP4OelQr

102◆2VB8wsVUoo 2017/09/04(月) 12:09:16.89ID:xP4OelQr

103◆2VB8wsVUoo 2017/09/04(月) 12:09:33.74ID:xP4OelQr

104◆2VB8wsVUoo 2017/09/04(月) 12:09:50.50ID:xP4OelQr

105◆2VB8wsVUoo 2017/09/04(月) 12:10:05.82ID:xP4OelQr

106◆2VB8wsVUoo 2017/09/04(月) 12:10:27.40ID:xP4OelQr

107◆2VB8wsVUoo 2017/09/04(月) 12:10:59.15ID:xP4OelQr

108◆2VB8wsVUoo 2017/09/04(月) 12:11:16.36ID:xP4OelQr

109132人目の素数さん2017/09/05(火) 06:09:36.65ID:PCfG056b
耳栓をしたら世界が変わってワロタ

110132人目の素数さん2017/09/19(火) 11:14:27.03ID:MBqwivnT
Reduce何て便利な関数が標準装備されていたんだな。

?Reduce
Reduce("+", c(1,2,3,4))
Reduce("+", c(1,2,3,4), accumulate = TRUE)
cumsum(c(1,2,3,4))
factorial(4)
Reduce('*',c(1,2,3,4))
Reduce("*", c(1,2,3,4), accumulate = TRUE)
cumprod(c(1,2,3,4))

Reduce("*", c(1,2,3,4), accumulate = TRUE,right=TRUE)
rev(Reduce("*", c(1,2,3,4), accumulate = TRUE,right=TRUE))
cumprod(rev(c(1,2,3,4)))

111132人目の素数さん2017/09/21(木) 19:43:24.15ID:KStrV4aS
96歳の女性をお看取り。肺も心臓も腎臓も機能低下していたが、年齢を考慮して死因は老衰にした。ふと96歳の平均余命はどれくらいかが気になって平成27年簡易生命表(女)を調べてみた。
http://www.mhlw.go.jp/toukei/saikin/hw/life/life15/dl/life15-07.pdfこういう表をみると自分で計算してみたくなる。諸関数の定義が厚労省の参考資料で解説されているのでRが使えれば計算は容易(ド底辺特殊シリツ医大卒には無理)

## 10万人あたりの死亡数から平均余命を算出
f=c(178,31,20,12,8,8,8,8,7,7,7,7,7,7,8,10,12,13,15,16,17,18,20,21,23,24,25,27,29,
30,31,32,34,37,39,41,43,46,50,56,62,68,73,79,86,94,103,113,123,134,145,159,174,
188,202,215,226,237,251,269,292,319,346,372,401,435,472,510,553,602,661,727,799,
870,952,1053,1180,1331,1503,1698,1914,2153,2419,2708,3011,3320,3628,3926,4227,
4513,4743,4896,4937,4867,4678,4385,4004,3557,3069,2567,2077,1623,1220,880,608,
960)
m=c(201,32,23,15,11,10,10,10,9,8,7,7,8,11,14,17,22,27,33,39,44,47,50,52,53,54,54,
54,54,56,58,59,62,65,70,73,76,79,85,93,103,113,122,131,145,159,176,195,215,237,
259,285,311,341,374,412,451,490,528,572,625,692,767,842,919,1002,1084,1166,1254,
1349,1452,1562,1668,1764,1875,2015,2182,2375,2586,2806,3038,3280,3506,3709,3887,
4024,4089,4088,4037,3944,3803,3555,3243,2887,2507,2122,1749,1403,1095,829,609,
434,298,198,127,180)

LE <-function(ndx,Y,N0=10^5){
n=length(ndx)
lx=numeric(n)
lx[1]=N0
for(i in 1:(n-1))
lx[i+1] <- lx[i] - ndx[i]
nqx=ndx/lx
nLx=numeric(n)
for(i in 1:n)
nLx[i] <- mean(c(lx[i],lx[i+1]))
nLx[n]=0
Tx=rev(cumsum(rev(nLx)))
le=Tx/lx
return(round(le[Y+1],1))
}

112132人目の素数さん2017/09/21(木) 19:43:57.06ID:KStrV4aS
>>356
実行すると
> LE(f,96)
[1] 3.4
96歳女性の平均余命は3年以上あるらしい。
死因は多臓器不全にする方が正しかったのだろうか?

> LE(m,100)
[1] 2.1

参考資料1 生命表諸関数の定義
http://www.mhlw.go.jp/toukei/saikin/hw/life/life15/dl/life15-08.pdf
に寿命中位数というのが定義されている。
半減期にあたる寿命中位数も計算してみた。グラフであたりをつけて内挿してもとめた。
手計算は面倒なのでRで線形回帰
x=c(89,90)
y=c(nLx[89+1],nLx[90+1]) # y = a + b*x
coef=lm(y~x)$coef
round((50000-coef[1])/coef[2],1) # 寿命中位数

> round((50000-coef[1])/coef[2],1) # 寿命中位数
89.3

参考資料中の解説

 我が国は、「平均寿命<寿命中位数」となっている。


が確認できた。

待機時間の多い職場だと、ふとした疑問をネットで調べて自分で勉強できて( ・∀・)イイ!!
ド底辺特殊シリツ医大卒のネット検索後の学習結果は>119を参照。

113132人目の素数さん2017/09/21(木) 20:07:38.28ID:g2Q4cSI1
どうでもいいけど、平成27年は完全生命表の年なのに、
なぜわざわざ簡易生命表を使うの?

114132人目の素数さん2017/09/23(土) 10:32:52.85ID:w2Jlt/XH
>>111
5〜6年ほど前にやってたな。

生命表と幾何ブラウン運動に基づく株価のモンテカルロシミュレーションで、今の貯蓄で死ぬまでにお金を使い切る確率を計算するの。

115132人目の素数さん2017/09/23(土) 19:46:43.88ID:lDVqqHRb
>>113
公表されていればurl希望

116132人目の素数さん2017/09/23(土) 21:22:11.76ID:Kb5kBoaF
>>115
第22回生命表(完全生命表)
ttp://www.mhlw.go.jp/toukei/saikin/hw/life/22th/index.html

117132人目の素数さん2017/09/23(土) 21:24:44.71ID:Kb5kBoaF
なお、最新版の平成28年の簡易生命表はこちら
ttp://www.mhlw.go.jp/toukei/saikin/hw/life/life16/index.html

118132人目の素数さん2017/09/24(日) 06:09:40.57ID:5AI1PnF/
>>117
>>116
ありがとうございました。

119132人目の素数さん2017/09/24(日) 11:13:15.02ID:o6BLfmLi
>>116
## 第22回生命表(完全生命表)
#http://www.mhlw.go.jp/toukei/saikin/hw/life/22th/index.html

F=c(178,32,20,12,8,8,8,8,7,7,7,7,7,7,8,10,12,13,15,16,17,19,20,22,23,24,25,27,28,30,31,32,34,36,39,41,42,45,50,56,
62,68,73,79,85,94,104,114,124,134,145,159,174,189,202,215,226,237,250,268,291,318,346,372,399,433,471,511,554,603,662,729,802,874,954,1053,1180,1332,1505,1699,1909,2143,2409,2701,3004,3310,3622,3938,4253,4531,4757,4918,5025,
5024,4876,4598,4132,3594,3025,2464,1941,1477,1085,769,525,345,217,132,76,42,22,11,5,2,1,0)

M=c(202,34,24,16,11,10,10,10,9,8,7,7,8,11,13,17,21,26,32,39,45,49,51,53,55,55,54,54,54,56,57,59,61,65,69,73,75,78,
84,93,103,113,122,131,144,159,176,195,215,236,257,283,310,340,373,411,450,488,525,568,620,688,764,839,910,994,
1081,1166,1256,1349,1450,1561,1675,1776,1885,2021,2185,2377,2594,2819,3046,3279,3504,3714,3900,4043,4116,4127,
4080,3973,3810,3580,3302,2967,2567,2123,1718,1352,1034,768,554,387,262,171,108,66,39,22,12,6,3,2,1)

LE <-function(ndx,Y,N0=10^5){
n=length(ndx)
lx=numeric(n)
lx[1]=N0
for(i in 1:(n-1))
lx[i+1] <- lx[i] - ndx[i]
nqx=ndx/lx
nLx=numeric(n)
for(i in 1:n)
nLx[i] <- mean(c(lx[i],lx[i+1]))
nLx[n]=0
Tx=rev(cumsum(rev(nLx)))
le=Tx/lx
return(round(le[Y+1],2))
}

120132人目の素数さん2017/09/24(日) 11:15:06.76ID:o6BLfmLi
>>117
最新版の平成28年の簡易生命表はこちら
f=c(198,29,19,12,9,7,7,6,5,6,6,7,7,7,8,9,11,12,13,14,16,19,22,24,25,25,25,25,25,26,28,29,31,34,37,41,45,48,50,54,59,66,72,78,84,92,101,112,123,136,149,162,175,186,196,207,220,237,256,275,294,313,334,360,
389,423,463,504,547,595,649,707,776,855,940,1043,1164,1308,1479,1667,1874,2102,2363,2651,2955,3265,3567,3874,4186,4484,4731,4908,5039,5067,4936,4635,4209,3697,3139,2574,2037,1554,1141,805,545,846)
m=c(194,31,21,14,10,9,9,8,7,7,7,7,8,10,13,17,21,26,31,38,44,48,50,51,51,51,53,53,54,56,57,58,60,63,66,70,74,79,83,89,97,106,116,127,140,155,171,190,211,233,255,278,302,328,359,395,434,478,525,572,622,
678,741,813,890,973,1062,1148,1233,1323,1419,1522,1633,1752,1874,2014,2176,2357,2553,2762,2985,3221,3459,3680,3875,4031,4123,4156,4123,4023,3874,3643,3349,3006,2629,2235,1843,1470,1131,837,593,401,258,157,89,90)
LE <-function(ndx,Y,N0=10^5){
n=length(ndx)
lx=numeric(n)
lx[1]=N0
for(i in 1:(n-1))
lx[i+1] <- lx[i] - ndx[i]
nqx=ndx/lx
nLx=numeric(n)
for(i in 1:n)
nLx[i] <- mean(c(lx[i],lx[i+1]))
nLx[n]=0
Tx=rev(cumsum(rev(nLx)))
le=Tx/lx
return(round(le[Y+1],1))
}
LE(f,96)
LE(m,100)

> LE(f,96) # 96歳女性の平均余命
[1] 3.3

> LE(m,100) #100歳男性の平均余命
[1] 1.9

微妙に結果が違うな。

121132人目の素数さん2017/09/27(水) 06:38:20.99ID:DzhNz5Oe
ド底辺特殊シリツ医大に
http://egg.2ch.net/test/read.cgi/hosp/1506113587/142
のようなスカトロ嗜好の医者がいるかを
1学年100人として無作為の順序で調査していき該当者が1人発見されたら調査終了とするという方法で調査する。
1例目が該当者であった。
この学年のスカトロ嗜好数の期待値を求めよ。

122132人目の素数さん2017/09/27(水) 06:39:04.84ID:DzhNz5Oe
lottery <- function(.N,.r,k=10^4){
f <-function(S,N=.N){
y=c(rep(1,S),rep(0,N-S))
Y=sample(y)
return(Y)
}
g <- function(y){
n=length(y)
for(i in 1:n){
if(!y[i]) i=i+1
else break
}
return(i)
}

xx=0:.N
SS=NULL
for(i in 1:k){
M=t(sapply(xx,f))
for(j in 1:.N){
if(g(M[j,])==.r) SS=c(SS,j-1)
}
}
print(quantile(SS,c(.025,.05,.50,.95,.975)))
c(mean=mean(SS))
}
lottery(100,1)

123132人目の素数さん2017/09/27(水) 06:39:43.99ID:DzhNz5Oe
> lottery(100,1)
2.5% 5% 50% 95% 97.5%
16 22 70 97 98
mean
66.31167

124132人目の素数さん2017/10/02(月) 22:03:04.10ID:jTVePLTX
lottery <- function(.N,.r,k=10^3){
f <-function(S,N=.N){
y=c(rep(1,S),rep(0,N-S))
Y=sample(y)
return(Y)
}
g <- function(y){
n=length(y)
for(i in 1:n){
if(!y[i]) i=i+1
else break
}
return(i)
}

xx=0:.N
SS=NULL
for(i in 1:k){
M=t(sapply(xx,f))
for(j in 1:(.N+1)){
if(g(M[j,])==.r) SS=c(SS,j-1)
}
}
print(quantile(SS,c(.025,.05,.50,.95,.975)))
c(mean=mean(SS))
}
lottery(100,1)

125132人目の素数さん2017/10/06(金) 08:15:11.45ID:zRAj2AHB
Linux版のRでもWindows版みたいに「入力:青字、出力:赤字」みたいなインターフェースにする設定ってある?

126132人目の素数さん2017/10/06(金) 19:21:19.13ID:nz6qypJh
>>125
coloroutパッケージ入れて、色を割り当てるとできるかも。
アーカイブに入っちゃっているけど、使える。
setOutputColors256()のヘルプを参照。

127132人目の素数さん2017/10/07(土) 05:15:49.03ID:8eYDySnN
version 3.4.2になってた。

128132人目の素数さん2017/10/07(土) 09:49:05.80ID:dws8gMbB
順番を表示させる、*1st,*2nd *3rd
> n=21
> paste(n,ifelse(0<n%%10 && n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th'),sep='')
[1] "21st"

129132人目の素数さん2017/10/07(土) 17:31:13.11ID:O5mMfwME
>>128
それだと、11が11stになっておかしいだろ

> library(toOrdinal)
> toOrdinal(21)
[1] "21st"
> sapply(1:20, toOrdinal)
[1] "1st" "2nd" "3rd" "4th" "5th" "6th" "7th" "8th" "9th" "10th"
[11] "11th" "12th" "13th" "14th" "15th" "16th" "17th" "18th" "19th" "20th"

130132人目の素数さん2017/10/07(土) 18:11:16.60ID:lKxjyQqB
11,12,13だけ除外かな

131132人目の素数さん2017/10/07(土) 22:26:53.02ID:W9vNhGzr
>>128
なんでそこで && を使うかなあ。&& の意味わかってる?
ifelse には & を使うことになるのが普通だし、そうでなければ if 〜 else のほうがよい。

132132人目の素数さん2017/10/08(日) 10:13:32.31ID:68bFYZIH
>>129
ご指摘とパッケージの御教示ありがとうございました。

133132人目の素数さん2017/10/08(日) 10:15:36.51ID:68bFYZIH
>>131
?&で勉強になりました。
御教示ありがとうございました。

134132人目の素数さん2017/10/08(日) 10:55:49.38ID:LDCLz3bM
>>130
toOrdinal のコードをみたら 11,12,13は別扱いしてあった。

if (tmp %in% c("11", "12", "13"))
tmp.suffix <- "th"

135132人目の素数さん2017/10/08(日) 11:36:14.74ID:LDCLz3bM
legendの中で使うべくone-linerにしたかったのだけど、ここでの指摘に従ってデバッグしたら

paste(n,ifelse(as.character(n%%100)%in%c('11','12','13'),'th',ifelse(0<n%%10 & n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th')),sep='')

こんなに長くなってしまった。

> f<-function(n)paste(n,ifelse(as.character(n%%100)%in%c('11','12','13'),'th',
+ ifelse(0<n%%10 & n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th')),sep='')
> prod(sapply(0:10000,f) == sapply(0:10000,toOrdinal::toOrdinal))
[1] 1

一応、パッケージtoOrdinalと0〜1万で合致を確認。

いろいろとご助言ありがとうございました。 <(_ _)>ペコリ

136132人目の素数さん2017/10/08(日) 11:55:19.96ID:68bFYZIH
パチンコのスレにこういう記載があった。

>バカは何故かパチンコで勝てると思い込んでいる
>バカは本来なら勝てるのに遠隔や不正で負けていると思い込んでいる
>バカは10分の1で10円が当たるクジより1万分の1で8000円が当たるクジの方が儲かると思ってる
>算数すらできないバカが必死に守って支えてきたのがパチンコ業界

これを読んでこんな問題を考えてみた。

宝くじAは1000本中100本が当たりで当たりは 100万円の賞金、

宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、

どちらも売り出し価格は同じなので計100本買うことにする。

Aを何本買うときに賞金の期待値が最大になるか、シミュレーションしてみよ。

137132人目の素数さん2017/10/08(日) 17:27:30.39ID:d4HM/FAz
分散が変わる?

138132人目の素数さん2017/10/08(日) 20:42:24.71ID:LDCLz3bM
シミュレーションすると

Na=1000 ; Nb=1000 # くじの本数
pa=0.1 ; pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
aa=0:n # 宝くじAを買う候補
A0=c(rep(1,Na*pa),rep(0,Na-Na*pa)) # 宝くじA 1:当たり 0:はずれ
B0=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb)) # 宝くじB

Get <- function(a){ #Aをa本買ったときのAの当たりとBのあたりでの賞金合計
sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb
}
k=10^3 #繰り返し回数
MAX=replicate(k,which.max(sapply(aa,Get))-1) # 賞金合計が最大になるAの購入数の配列
hist(MAX)

期待値がA,BでかわらないからMAXの分布は一様分布になると思ったのだが、

こんな分布になってしまって自分でも混乱している。

http://i.imgur.com/Yo55mvc.png

お知恵を拝借したい。

139132人目の素数さん2017/10/08(日) 22:12:35.95ID:COe73DTR
よくわからん。rbinomとか使えば?

140132人目の素数さん2017/10/08(日) 22:47:20.08ID:7IiqlruN
>>139
俺もそう思った

141132人目の素数さん2017/10/08(日) 23:23:18.27ID:7IiqlruN
自分なら、こんな感じ
> a <- function(i){sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)}
> b <- function(){which.max(sapply(1:100, a))}
1000回のシミュレーションだと
> res <- replicate(1000, b())
> hist(res)
同じ結果だよ

142132人目の素数さん2017/10/08(日) 23:29:05.07ID:7IiqlruN
A, Bを少なくとも1本含む100本じゃなかったので、
> b <- function(){which.max(sapply(1:100, a))}

> b <- function(){which.max(sapply(0:100, a))}
だな。すまん
> a <- function(i){sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)}
> b <- function(){which.max(sapply(0:100, a))}
> res <- replicate(1000, b())
> hist(res,breaks = 20)
結果は同じ

143132人目の素数さん2017/10/09(月) 00:19:33.56ID:1+Fkporl
>139-142

簡潔なコードありがとうございます。

AもBも賞金の期待値がどちらも10万なので
Aを何本買おうが合計の期待値は同じかなと思いつつ、sample関数を使ってシミュレーションしたのだけれど
一様分布にはならないのでシミュレーションの間違い、
もしくは自分が理解していないsample関数の特性を検出しただけなのか、と思っておりました。
rbinomでのシミュレーション結果とヒストグラムが一致して、一様分布になるという思い込みが間違いだと確信できました。

勉強になりました。ありがとうございました。

144132人目の素数さん2017/10/09(月) 00:27:22.71ID:1+Fkporl
>>142
大したことではないけど

0から始まるので本数とindexを補正して

b <- function(){which.max(sapply(0:100, a))-1}

の方が正確だと思います。

145132人目の素数さん2017/10/09(月) 01:00:19.99ID:iikzqEdl
合計の期待値の標準偏差はAの本数が増えると減るね

146132人目の素数さん2017/10/09(月) 01:05:40.45ID:1+Fkporl
宝くじをA,B併せて100本買ったときの賞金の期待値は

Award <- function(i){ #i:Aを買った本数
sum(dbinom(0:i,i,0.1)*(0:i)*100)+sum(dbinom(0:(100-i),100-i,0.01)*(0:(100-i))*1000)
}
sapply(0:100,Award)

> sapply(0:100,Award)
[1] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
....
[86] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000

と買い方に無関係にも思えるんだなぁ。

147132人目の素数さん2017/10/09(月) 02:00:05.31ID:1+Fkporl
>>145
>合計の期待値の標準偏差はAの本数が増えると減る
ことを検証してみた。

a <- function(i) sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)
plot(0:100,sapply(0:100,function(x) sd(replicate(10^3,a(x)))))

http://i.imgur.com/oPFtMXe.png

148132人目の素数さん2017/10/09(月) 06:15:43.03ID:1+Fkporl
別の話題提供。

 ゴルゴ13は100発100中

 ゴルゴ14は10発10中

 ゴルゴ15は1発1中

とする。

各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか?

149132人目の素数さん2017/10/09(月) 15:22:52.36ID:di3eoSv3
どうでもいいんだが、
which.max()は、同値の最大が複数あったとき、最初のインデックスを戻すので、
2番目以降の最大が無視されていることについて、その取り扱いを考慮した方が良いと思うぞ。
> which.max(c(0,1,0,0,0,0,0,0,1))
[1] 2

150132人目の素数さん2017/10/09(月) 17:42:12.61ID:1+Fkporl
>>149
ご指摘ありがとうございます。

最大値を与えるAの購入数が複数ある場合を考えてスクリプトを改変してみました。

a <- function(i){sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)}
which.max2 <- function(x){which(x==max(x))-1}
b <- function(){
b1=NULL
b1=c(b1,which.max2(sapply(0:100, a)))
return(b1)
}
b2=NULL
for(j in 1:1000){
b2=c(b2,b())
}
hist(b2,breaks=20)

151132人目の素数さん2017/10/09(月) 19:43:31.10ID:1+Fkporl
length(b2)

とすると1000を超えているから同値の最大が複数あったことがわかりました。

ありがとうございました。

このスレは勉強になるなぁ。

152132人目の素数さん2017/10/09(月) 22:05:49.50ID:1+Fkporl
どうでもいいことだけど、

hist(b2,breaks=20,col=sample(colors(),2))

とすると毎回配色の変わるツートンカラーのヒストグラムがでて面白い。

153132人目の素数さん2017/10/10(火) 18:01:18.51ID:S7i/7cX1
>>126
125です。レス遅くなったけど、やりたいことができました。ありがとう。

154132人目の素数さん2017/10/11(水) 20:55:10.55ID:DiiGwH6/
array(1:48, dim=c(6, 4, 2))
を,3行2列おきに平均して
,,1
[,1] [,2]
[1,] 5 17
[2,] 8 20
,,2
[,1] [,2]
[1,] 29 41
[2,] 32 44
としたいんですが,効率的な方法はないでしょうか。
実際のデータは,10000x10000x400なんです。

155132人目の素数さん2017/10/11(水) 22:16:13.22ID:bFq2a2U4
>>154
>3行2列
さっぱりわからん。
例えば、結果の[1,1,1]の5は、どれとどれの平均?

156132人目の素数さん2017/10/11(水) 22:25:55.39ID:HJHHXjdz
>>155
6行4列の左上の6つ1,2,3,7,8,9の平均5だと思う。

157132人目の素数さん2017/10/11(水) 22:59:20.81ID:HJHHXjdz
>>154
効率的かどうかはわからんが、やってみた

(d=array(1:48, dim=c(6, 4, 2)))

a=3
b=2
l=6
m=4
n=2

re=NULL
for(k in 1:n){
for(j in 1:(m/b)){
for(i in 1:(l/a)){
re=append(re,mean(d[1:3+3*(i-1),1:2+2*(j-1),k]))
}
}
}
array(re,dim=c(m/b,l/a,n))

158132人目の素数さん2017/10/11(水) 23:00:23.97ID:UQm0Ph1Q
>>155
>>154です。
説明が悪くてすみません。
>>156の通りです。

159132人目の素数さん2017/10/11(水) 23:02:02.80ID:UQm0Ph1Q
>>157
ありがとうございます。
手元にデータがないので明日試してみます。

160132人目の素数さん2017/10/11(水) 23:02:32.81ID:bFq2a2U4
なるほど、じゃあその5,8,17,20はこういうことか。
> a <- array(1:48, dim=c(6, 4, 2))
> b <- expand.grid(1:(dim(a)[1] %/% 3), 1:(dim(a)[2] %/% 2))
> f <- function(i, j, k=1){mean(a[(3*i-2):(3*i),(2*j-1):(2*j),k])}
> mapply(f, b[,1], b[,2])
[1] 5 8 17 20

161132人目の素数さん2017/10/11(水) 23:04:06.69ID:bFq2a2U4
あぁ、すでに >>157 が回答を書いていたorz

162132人目の素数さん2017/10/11(水) 23:06:59.19ID:HJHHXjdz
> d=array(1:48, dim=c(6, 4, 2))

> a=3
> b=2
> l=6
> m=4
> n=2
> re=NULL
> for(k in 1:n){
+ for(j in 1:(m/b)){
+ for(i in 1:(l/a)){
+ re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),k]))
+ }
+ }
+ }
> array(re,dim=c(m/b,l/a,n))
, , 1

[,1] [,2]
[1,] 5 17
[2,] 8 20

, , 2

[,1] [,2]
[1,] 29 41
[2,] 32 44

データは10000x10000x400なら
l=10000
m=10000
n=400
でよくね?

163132人目の素数さん2017/10/12(木) 07:10:38.11ID:2sFs7Gsb
最初の1行で

> d=array(1:(10000*10000*400))
Error: cannot allocate vector of size 298.0 Gb

エラーがでたw

164132人目の素数さん2017/10/12(木) 07:12:09.84ID:2sFs7Gsb
正しいエラーwはこっち

> d=array(1:(10000*10000*400),dim=c(6, 4, 2))
Error: cannot allocate vector of size 298.0 Gb

165132人目の素数さん2017/10/12(木) 07:21:02.01ID:2sFs7Gsb
そもそも、10000x10000x400だと10000行は3行にならないな。最終行を無視して9999行にするのかな?

166132人目の素数さん2017/10/12(木) 14:48:46.91ID:sG1WSodN
>160を参考に>157のforだらけのスパゲティー・ソースを(ちょっと一般化して)改変しました。
(変数名の重複を避けるため>160の変数名は変更してあります。)

> a=3
> b=2
> c=1
>
> l=6
> m=4
> n=2
>
> A <- array(1:48, dim=c(l,m,n))
> B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c))
> f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])}
> re=mapply(f, B[,1], B[,2],B[,3])
> array(re,dim=c(m%/%b,l%/%a,n%/%c))
, , 1

[,1] [,2]
[1,] 5 17
[2,] 8 20

, , 2

[,1] [,2]
[1,] 29 41
[2,] 32 44

>

上級者のコードを読むのは勉強になります。(深謝)

167132人目の素数さん2017/10/12(木) 16:33:51.03ID:f/ETd6de
先輩方,>>154です。
巨細にわたってご教示・ご議論いただきありがとうございました。
>>166のコードを参考にして,無事処理することができました。
実行速度も申し分なく,効率的に作業が進められそうです。
今回は私には目から鱗な解法で,勉強になりました。
質問して良かったと思っています。。本当にありがとうございました。

10,000行の扱いですが,3行ごとに処理して1行無視してもよいし,
4から7行程度で処理してもよく,フレキシブルです。
今回は4行にして余りなしにしました。

168132人目の素数さん2017/10/12(木) 18:43:54.51ID:2sFs7Gsb
処理速度を比べてみた。

fs<-function(){ #スパゲッティ・ソース
d=array(1:(l*m*n),dim=c(l,m,n))
re=NULL
for(k in 1:(n/c)){
for(j in 1:(m/b)){
for(i in 1:(l/a)){
re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),1:c+c*(k-1)]))
}}}
array(re,dim=c(m/b,l/a,n))
}

fe <- function(){ # Expertソース
A <- array(1:(l*m*n), dim=c(l,m,n))
B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c))
f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])}
re=mapply(f, B[,1], B[,2],B[,3])
array(re,dim=c(m%/%b,l%/%a,n%/%c))
}

a=4
b=2
c=1
l=400
m=200
n=2

> system.time(fs())
user system elapsed
2.17 0.00 2.18
> system.time(fe())
user system elapsed
1.14 0.00 1.14

>160氏のスキルに改めて感服。

169132人目の素数さん2017/10/14(土) 23:02:55.66ID:WIgy7a2M
# 既約剰余類群 (Z/nZ)*
kjrg <- function(n){
nn=0:(n-1)
f=function(x,y) (x*y)%%n
names(nn)=paste0('n',nn)
z=outer(nn,nn,f)
idx=which(gmp::gcd(1:(n-1),n)==1)+1
return(z[idx,idx])
}
kjrg(10)
kjrg(15)
kjrg(24)
> kjrg(24)
n1 n5 n7 n11 n13 n17 n19 n23
n1 1 5 7 11 13 17 19 23
n5 5 1 11 7 17 13 23 19
n7 7 11 1 5 19 23 13 17
n11 11 7 5 1 23 19 17 13
n13 13 17 19 23 1 5 7 11
n17 17 13 23 19 5 1 11 7
n19 19 23 13 17 7 11 1 5
n23 23 19 17 13 11 7 5 1

170132人目の素数さん2017/10/15(日) 13:21:32.63ID:GlVyhxtn
# 既約剰余類群での加減乗除
n=24
n_star=which(gmp::gcd(1:(n-1),n)==1) ; n_star
tasu <- function(x,y) (x+y)%%n
hiku <- function(x,y) (x-y)%%n # row - col
kake <- function(x,y) (x*y)%%n
g=function(x) n_star[which(x==1)]
.M=outer(n_star,n_star,kake)
G=apply(.M,2,g)
gyaku <- function(x) n_star[which(G==(x%%n))]
waru <- function(x,y) (x*gyaku(y))%%n # row ÷ col

xx=yy=c(0,n_star)
names(xx)=paste0('x',c(0,n_star))
names(yy)=paste0('y',c(0,n_star))
outer(xx,yy,tasu) # x + y
outer(xx,yy,hiku) # x - y
outer(xx,yy,kake) # x * y
X=Y=n_star #outer(X,Y,waru) # WRONG!!
a=expand.grid(X,Y)
b=matrix(mapply(waru,a[,1],a[,2]),length(X))
rownames(b)=paste0('x',n_star)
colnames(b)=paste0('y',n_star)
b # x / y

171132人目の素数さん2017/10/15(日) 14:03:02.02ID:GlVyhxtn
名前も一致してidenticalなんだな。
> x=1:10
> y=1:10
> identical(x,y)
[1] TRUE
> names(x)=LETTERS[1:10]
> identical(x,y)
[1] FALSE
> names(y)=LETTERS[1:10]
> identical(x,y)
[1] TRUE
> names(y)=letters[1:10]
> identical(x,y)
[1] FALSE

allの方は内容の全体が一致していれば順番は入れ替わってもTRUEなんだな。
> x=1:10
> y=1:10
> all(x,y)
[1] TRUE
> names(x)=LETTERS[1:10]
> all(x,y)
[1] TRUE
> names(y)=letters[1:10]
> all(x,y)
[1] TRUE
> all(x,sample(y))
[1] TRUE

172132人目の素数さん2017/10/16(月) 01:11:00.56ID:NsB0m+v3
beki <- function(x,n,p){ # x^n%%p
if(n==0) return(1)
if(n==1) return(x%%p)
re=numeric(n)
re[1]=x
for(i in 1:(n-1)) re[i+1] = (x*re[i])%%p
return(re[n])
}

p=1009 # フェルマーの小定理
xx=0:(p-1)
n=p-1
y=sapply(xx,function(x)beki(x,n,p))
summary(y)

173132人目の素数さん2017/10/16(月) 07:46:39.05ID:NsB0m+v3
90歳女性のお看取り

>120のコードから女性の平均余命をグラフにしてみた
Age=85:105
(Life_Expectancy=sapply(Age,function(x) LE(f,x)))
plot(Age,Life_Expectancy,pch=19,col=sample(colors()))
abline(h=5,col='gray',lty=2)

http://i.imgur.com/38x50fw.png

平均余命が5年を切るのは92歳なので、老衰にせずに
主治医指定の病名で診断書を作成。
カルテ記載のしっかりした患者の診療は負担なくこなせて(・∀・)イイ!!

174132人目の素数さん2017/10/16(月) 22:16:44.56ID:ACageh59
outerって使えない場面があるよね?
剰余系での加減乗除での実際
> n=5 # prime number
> nn=1:(n-1)
> tasu <- function(x,y) (x+y)%%n
> hiku <- function(x,y) (x-y)%%n # row - col
> kake <- function(x,y) (x*y)%%n
> g=function(x) nn[which(x==1)]
> .M=outer(nn,nn,kake)
> G=apply(.M,2,g)
> gyaku <- function(x) nn[which(G==(x%%n))]
> waru <- function(x,y) (x*gyaku(y))%%n # row / col
> waru(3,2)
[1] 4
> xx=yy=c(0,nn)
> names(xx)=paste0('x',c(0,nn))
> names(yy)=paste0('y',c(0,nn))
> outer(xx,yy,tasu) # x + y
y0 y1 y2 y3 y4
x0 0 1 2 3 4
x1 1 2 3 4 0
x2 2 3 4 0 1
x3 3 4 0 1 2
x4 4 0 1 2 3
> outer(xx,yy,hiku) # x - y
y0 y1 y2 y3 y4
x0 0 4 3 2 1
x1 1 0 4 3 2
x2 2 1 0 4 3
x3 3 2 1 0 4
x4 4 3 2 1 0
> outer(xx,yy,kake) # x * y
y0 y1 y2 y3 y4
x0 0 0 0 0 0
x1 0 1 2 3 4
x2 0 2 4 1 3
x3 0 3 1 4 2
x4 0 4 3 2 1
> X=Y=nn
> outer(X,Y,waru) # WRONG!!
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] NA NA NA NA
[3,] NA NA NA NA
[4,] NA NA NA NA
> a=expand.grid(X,Y)
> b=matrix(mapply(waru,a[,1],a[,2]),length(X))
> rownames(b)=paste0('x',nn)
> colnames(b)=paste0('y',nn)
> b # x / y
y1 y2 y3 y4
x1 1 3 2 4
x2 2 1 4 3
x3 3 4 1 2
x4 4 2 3 1

175132人目の素数さん2017/10/16(月) 23:28:56.68ID:FS7yMelZ
>>174
そりゃあ、outer に渡す関数がまずいんだろ。
outer に渡す関数は、完全にベクトル対応になっていなければならないんじゃないの?
でなければ、わざわざ outer 使うこともなかろう。

176132人目の素数さん2017/10/17(火) 06:58:05.54ID:ps5i5Wv8
>>175
レスありがとうございます。
自作関数から自作関数を呼び出しているのが原因かと思っていましたが、
ベクトル対応していないとはこういうことでしょうか?

加減乗は引数のどちらをベクトルにしても動作
> kake(1:4,1)
[1] 1 2 3 4
> kake(1,1:4)
[1] 1 2 3 4
> tasu(1:4,1)
[1] 2 3 4 0
> tasu(1,1:4)
[1] 2 3 4 0
> hiku(1:4,1)
[1] 0 1 2 3
> hiku(1,1:4)
[1] 0 4 3 2

除算は第2引数をベクトルにすると誤動作

> waru(1:4,1)
[1] 1 2 3 4
> waru(4,1:4)
[1] 4 1

177132人目の素数さん2017/10/17(火) 22:43:30.39ID:yEir4ohi
>>176
outer から呼び出されたときの引数を見てみればわかるよ。たぶん1度しか呼ばれず、できあがりの行列の要素すべてを計算するための引数がベクトルで渡されているはずだから。

178名無しさん@そうだ選挙に行こう! Go to vote!2017/10/22(日) 07:51:52.36ID:/ES4CbkI
>>177
ありがとうございます。

outerのソースをみたら、こんな記述があったので、X,Yを同時にベクトルで与えても動作する関数でないとダメみたいでした。

FUN <- match.fun(FUN)
Y <- rep(Y, rep.int(length(X), length(Y)))
if (length(X))
X <- rep(X, times = ceiling(length(Y)/length(X)))
FUN(X, Y, ...)

179名無しさん@そうだ選挙に行こう! Go to vote!2017/10/22(日) 08:28:19.35ID:/ES4CbkI
# greatest common divisor
gcd <- function(x,y) ifelse(x|y,ifelse(x*y,max(intersect((1:abs(x))[abs(x)%%(1:abs(x))==0],(1:abs(y))[abs(y)%%(1:abs(y))==0])),max(abs(x),abs(y))),NA)

非正整数対応すると長くなるなぁ

> gcd(24,15)
[1] 3
> gcd(16,15)
[1] 1
> gcd(5,0)
[1] 5
> gcd(-24,15)
[1] 3
> gcd(-24,-15)
[1] 3
> gcd(0,0)
[1] NA

180名無しさん@そうだ選挙に行こう! Go to vote!2017/10/22(日) 09:17:13.05ID:/ES4CbkI
ユークリッドの互除法をつかうとこんな感じ

> GCD <- function(x,y){
+ if(round(x)!=x | round(y)!=y) stop('Not integer')
+ if(!x&!y) return(NA)
+ a=max(c(abs(x),abs(y)))
+ b=min(c(abs(x),abs(y)))
+ if(!a*b) return(a)
+ r=integer()
+ r[1]=a
+ r[2]=b
+ i=1
+ r[i+2]=r[i]%%r[i+1]
+ while(r[i+2]!=0 & r[i+2]!=1){
+ i=i+1
+ r[i+2]=r[i]%%r[i+1]
+ }
+ return(ifelse(r[i+2]==0,r[i+1],1))
+ }
> GCD(0,0)
[1] NA
> GCD(24,15)
[1] 3
> GCD(16,15)
[1] 1
> GCD(5,0)
[1] 5
> GCD(24,-15)
[1] 3
> GCD(-24,-15)
[1] 3
> GCD(2.3,4)
Error in GCD(2.3, 4) : Not integer

181名無しさん@そうだ選挙に行こう! Go to vote!2017/10/22(日) 10:07:52.18ID:tFigTWxc
文系研究者のワイ、困惑
ディープな世界やでほんま

182132人目の素数さん2017/10/25(水) 20:38:52.07ID:A26rYIkF
#.n発r中の狙撃手が.N発狙撃するときの命中数を返す
Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
print(summary(SS))
invisible(SS)
}

183◆2VB8wsVUoo 2017/10/26(木) 04:11:35.13ID:jT09z118

184◆2VB8wsVUoo 2017/10/26(木) 04:11:58.98ID:jT09z118

185◆2VB8wsVUoo 2017/10/26(木) 04:12:18.59ID:jT09z118

186◆2VB8wsVUoo 2017/10/26(木) 04:12:35.94ID:jT09z118

187◆2VB8wsVUoo 2017/10/26(木) 04:12:53.44ID:jT09z118

188◆2VB8wsVUoo 2017/10/26(木) 04:13:10.90ID:jT09z118

189◆2VB8wsVUoo 2017/10/26(木) 04:13:28.27ID:jT09z118

190◆2VB8wsVUoo 2017/10/26(木) 04:13:46.40ID:jT09z118

191◆2VB8wsVUoo 2017/10/26(木) 04:14:03.59ID:jT09z118

192◆2VB8wsVUoo 2017/10/26(木) 04:14:22.57ID:jT09z118

193132人目の素数さん2017/10/26(木) 18:17:25.90ID:tHqsRRzo
こういう問題を考えるのは楽しいね。

ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか?

ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか?

194132人目の素数さん2017/11/01(水) 12:40:01.83ID:UKV+bz4p
ベクトルは再利用されるんだな


> c(1,2,1,2)==c(1,2)
[1] TRUE TRUE TRUE TRUE
> all(c(1,2,1,2)==c(1,2))
[1] TRUE
> equal <- function(x,y) sum((x-y)^2)==0 & length(x)==length(y)
> equal(c(1,2,1,2),c(1,2))
[1] FALSE

195132人目の素数さん2017/11/03(金) 08:27:39.47ID:7YwYHhPi
> x= list(length=5)
> x
$length
[1] 5

> y=list()
> length(y)=4
> y
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

> z=list(5)
> z
[[1]]
[1] 5

> w = vector(mode="list" , length= 5)
> w
[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

[[5]]
NULL

196132人目の素数さん2017/11/09(木) 08:41:04.22ID:SxG1YX2j
Brunner-Munzel検定って比較する両群に重なりがないとp値がNAになるんだが、
これってバグなのか?



> x=1:5
> y=6:10
> lawstat::brunner.munzel.test(x,y)

Brunner-Munzel Test

data: x and y
Brunner-Munzel Test Statistic = Inf, df = NaN, p-value = NA
95 percent confidence interval:
NaN NaN
sample estimates:
P(X<Y)+.5*P(X=Y)
1

> x=1:6
> y=6:11
> lawstat::brunner.munzel.test(x,y)

Brunner-Munzel Test

data: x and y
Brunner-Munzel Test Statistic = 24.749, df = 10, p-value =
2.651e-10
95 percent confidence interval:
0.9423463 1.0298759
sample estimates:
P(X<Y)+.5*P(X=Y)
0.9861111

197132人目の素数さん2017/11/14(火) 22:02:36.50ID:Jlt4ZA5T
突然ですが質問があります。最近Rソフトを利用し始めた者のです。単刀直入に、
ボロノイ図を制作した場合に分割された面積の平均を出力するには、どうすればよいのでしょうか?

198132人目の素数さん2017/11/14(火) 22:06:32.42ID:Jlt4ZA5T
現状としては住所をグーグルマップの緯度経度に変換し、分割まではできています。
宜しければどなたかお力添えをいただけないでしょうか?

199◆2VB8wsVUoo 2017/11/15(水) 00:58:17.28ID:WZuPK5Ir

200◆2VB8wsVUoo 2017/11/15(水) 00:58:38.82ID:WZuPK5Ir

201◆2VB8wsVUoo 2017/11/15(水) 00:58:56.67ID:WZuPK5Ir

202◆2VB8wsVUoo 2017/11/15(水) 00:59:14.81ID:WZuPK5Ir

203◆2VB8wsVUoo 2017/11/15(水) 00:59:32.78ID:WZuPK5Ir

204◆2VB8wsVUoo 2017/11/15(水) 00:59:50.57ID:WZuPK5Ir

205◆2VB8wsVUoo 2017/11/15(水) 01:00:09.89ID:WZuPK5Ir

206◆2VB8wsVUoo 2017/11/15(水) 01:00:29.29ID:WZuPK5Ir

207◆2VB8wsVUoo 2017/11/15(水) 01:00:47.81ID:WZuPK5Ir

208◆2VB8wsVUoo 2017/11/15(水) 01:01:09.56ID:WZuPK5Ir

209132人目の素数さん2017/11/15(水) 17:11:14.46ID:eZCzfwkb
>>198
SpatialPolygonsクラスまでできるといると仮定して、
各ポリゴンの中にareaスロットがあるから、
それを抽出して、その平均を取れば良い。
ボロノイ図がSpatialPolygonsクラスになっていないなら、
まずは変換。

210◆2VB8wsVUoo 2017/11/15(水) 20:30:49.23ID:WZuPK5Ir

211◆2VB8wsVUoo 2017/11/15(水) 20:31:06.31ID:WZuPK5Ir

212◆2VB8wsVUoo 2017/11/15(水) 20:31:21.95ID:WZuPK5Ir

213◆2VB8wsVUoo 2017/11/15(水) 20:31:40.54ID:WZuPK5Ir

214◆2VB8wsVUoo 2017/11/15(水) 20:31:58.15ID:WZuPK5Ir

215◆2VB8wsVUoo 2017/11/15(水) 20:32:16.33ID:WZuPK5Ir

216◆2VB8wsVUoo 2017/11/15(水) 20:32:32.65ID:WZuPK5Ir

217◆2VB8wsVUoo 2017/11/15(水) 20:32:48.32ID:WZuPK5Ir

218◆2VB8wsVUoo 2017/11/15(水) 20:33:05.48ID:WZuPK5Ir

219◆2VB8wsVUoo 2017/11/15(水) 20:33:21.41ID:WZuPK5Ir

220132人目の素数さん2017/11/15(水) 21:12:29.48ID:9zgzmHZc
ボロノイ図の各領域の面積の平均値って、単に
描画領域の総面積÷点数 じゃないの?
何か暗黙の前提があるの?

221◆2VB8wsVUoo 2017/11/16(木) 06:01:20.12ID:6ldUKvsQ

222◆2VB8wsVUoo 2017/11/16(木) 06:01:37.22ID:6ldUKvsQ

223◆2VB8wsVUoo 2017/11/16(木) 06:01:56.08ID:6ldUKvsQ

224◆2VB8wsVUoo 2017/11/16(木) 06:02:14.27ID:6ldUKvsQ

225◆2VB8wsVUoo 2017/11/16(木) 06:02:31.67ID:6ldUKvsQ

226◆2VB8wsVUoo 2017/11/16(木) 06:02:48.95ID:6ldUKvsQ

227◆2VB8wsVUoo 2017/11/16(木) 06:03:06.96ID:6ldUKvsQ

228◆2VB8wsVUoo 2017/11/16(木) 06:03:24.98ID:6ldUKvsQ

229◆2VB8wsVUoo 2017/11/16(木) 06:03:42.58ID:6ldUKvsQ

230◆2VB8wsVUoo 2017/11/16(木) 06:04:00.50ID:6ldUKvsQ

231132人目の素数さん2017/11/17(金) 20:18:17.65ID:5lHMDqma
Stan って 離散量をparameter 扱いできないので

m ? categorical(θ) としても m の分布を出せない。

JAGSだと

m 〜 dcat(θ)

とすると mの分布がだせる。

この理解でいいかな?

232◆2VB8wsVUoo 2017/11/17(金) 20:41:11.51ID:Y8c01QBt

233◆2VB8wsVUoo 2017/11/17(金) 20:41:27.83ID:Y8c01QBt

234◆2VB8wsVUoo 2017/11/17(金) 20:41:43.62ID:Y8c01QBt

235◆2VB8wsVUoo 2017/11/17(金) 20:41:58.17ID:Y8c01QBt

236◆2VB8wsVUoo 2017/11/17(金) 20:42:13.43ID:Y8c01QBt

237◆2VB8wsVUoo 2017/11/17(金) 20:42:28.14ID:Y8c01QBt

238◆2VB8wsVUoo 2017/11/17(金) 20:42:43.14ID:Y8c01QBt

239◆2VB8wsVUoo 2017/11/17(金) 20:42:59.36ID:Y8c01QBt

240◆2VB8wsVUoo 2017/11/17(金) 20:43:14.19ID:Y8c01QBt

241◆2VB8wsVUoo 2017/11/17(金) 20:43:30.82ID:Y8c01QBt

242132人目の素数さん2017/11/20(月) 17:47:00.70ID:ptVXRSQM
以前、ボロノイ図について質問させていただいた初心者です。
質問後色々と試行錯誤したのですがdir.areaを出せず悩んでおります。
library(deldir)
> library(ggmap)
> df <- read.csv("C:/Users/***/Desktop/***/***.csv",header=T)
> loc <- c(min(df$lon),min(df$lat),max(df$lon),max(df$lat))
> GM <-ggmap(get_map(location=loc,zoom=11,source="google"))
> GP <- geom_point(data=df,aes(x=lon,y=lat,colour=factor(type)),size=3)
> dd <- deldir(df$lon,df$lat)
>GS <-geom_segment(data=dd$dirsgs,aes(x=x1,y=y1,xend=x2,yend=y2),size=0.5,linetype=1)
> GM+GP+GS
と入力しボロノイ図は出せているのですが、

243132人目の素数さん2017/11/20(月) 17:52:19.12ID:ptVXRSQM
summaryと入力しても
function (object, ...)
UseMethod("summary")
<bytecode: 0x06b5eef0>
<environment: namespace:base>
と表示されます。どの様に入力または再入力すればdir.areaが表示させられるでしょうか?
私自身初めてで、全く分からない状態でRを使用しているので、宜しければ分かり易く教えて頂けると幸いです。

244◆2VB8wsVUoo 2017/11/20(月) 21:05:10.85ID:yjl62pX+

245◆2VB8wsVUoo 2017/11/20(月) 21:05:28.54ID:yjl62pX+

246◆2VB8wsVUoo 2017/11/20(月) 21:05:44.92ID:yjl62pX+

247◆2VB8wsVUoo 2017/11/20(月) 21:06:01.26ID:yjl62pX+

248◆2VB8wsVUoo 2017/11/20(月) 21:06:18.27ID:yjl62pX+

249◆2VB8wsVUoo 2017/11/20(月) 21:06:35.26ID:yjl62pX+

250◆2VB8wsVUoo 2017/11/20(月) 21:06:52.55ID:yjl62pX+

251◆2VB8wsVUoo 2017/11/20(月) 21:07:09.97ID:yjl62pX+

252◆2VB8wsVUoo 2017/11/20(月) 21:07:27.71ID:yjl62pX+

253◆2VB8wsVUoo 2017/11/20(月) 21:07:44.86ID:yjl62pX+

254132人目の素数さん2017/11/22(水) 10:59:21.45ID:1eSvd7ei
>>242
検証用のデータがないと誰も追試しないと思う。
>154は模擬データがあったのでやる気になった。

255◆2VB8wsVUoo 2017/11/22(水) 12:06:07.41ID:j3iP9uSb

256◆2VB8wsVUoo 2017/11/22(水) 12:06:26.73ID:j3iP9uSb

257◆2VB8wsVUoo 2017/11/22(水) 12:06:45.80ID:j3iP9uSb

258◆2VB8wsVUoo 2017/11/22(水) 12:07:02.42ID:j3iP9uSb

259◆2VB8wsVUoo 2017/11/22(水) 12:07:21.82ID:j3iP9uSb

260◆2VB8wsVUoo 2017/11/22(水) 12:07:41.50ID:j3iP9uSb

261◆2VB8wsVUoo 2017/11/22(水) 12:08:00.58ID:j3iP9uSb

262◆2VB8wsVUoo 2017/11/22(水) 12:08:22.66ID:j3iP9uSb

263◆2VB8wsVUoo 2017/11/22(水) 12:08:39.11ID:j3iP9uSb

264◆2VB8wsVUoo 2017/11/22(水) 12:08:57.12ID:j3iP9uSb

265132人目の素数さん2017/11/24(金) 16:46:01.01ID:S4G8MO/P
関数のオプションの取り得る値を全て表示する方法って知らない?
例えばcor関数のmethodオプションの取り得る値は"pearson", "kendall", "spearman"みたいに。
そういうのは無いのかな?
初めて使う関数のオプションいちいちググるの面倒くさい。

266◆2VB8wsVUoo 2017/11/24(金) 18:45:27.78ID:7RwNGOaZ

267◆2VB8wsVUoo 2017/11/24(金) 18:45:42.78ID:7RwNGOaZ

268◆2VB8wsVUoo 2017/11/24(金) 18:45:58.90ID:7RwNGOaZ

269◆2VB8wsVUoo 2017/11/24(金) 18:46:16.24ID:7RwNGOaZ

270◆2VB8wsVUoo 2017/11/24(金) 18:46:33.09ID:7RwNGOaZ

271◆2VB8wsVUoo 2017/11/24(金) 18:46:49.63ID:7RwNGOaZ

272◆2VB8wsVUoo 2017/11/24(金) 18:47:03.93ID:7RwNGOaZ

273◆2VB8wsVUoo 2017/11/24(金) 18:47:38.91ID:7RwNGOaZ

274◆2VB8wsVUoo 2017/11/24(金) 18:47:57.22ID:7RwNGOaZ

275◆2VB8wsVUoo 2017/11/24(金) 18:48:16.46ID:7RwNGOaZ

276132人目の素数さん2017/11/24(金) 19:33:11.50ID:aKG4T9LS
>>265
ヘルプ表示じゃ駄目なのかい?

277132人目の素数さん2017/11/24(金) 20:02:53.10ID:YIwuPhoP
>>265
> args(cor)
function (x, y = NULL, use = "everything", method = c("pearson",
"kendall", "spearman"))
NULL

278132人目の素数さん2017/11/24(金) 22:49:54.91ID:8tVkv6KQ
>>276
前にヘルプで見ようと思ったら書いてなかった覚えが…。マニアックな関数だったんで、そんなもんかな〜と諦めてた。やっぱ基本はヘルプか。
>>277
argsは見れる関数と見れない関数があるのであんまり役に立たない印象。

279◆2VB8wsVUoo 2017/11/25(土) 00:22:59.73ID:kb2pRtxG

280◆2VB8wsVUoo 2017/11/25(土) 00:23:21.36ID:kb2pRtxG

281◆2VB8wsVUoo 2017/11/25(土) 00:23:46.66ID:kb2pRtxG

282◆2VB8wsVUoo 2017/11/25(土) 00:24:03.98ID:kb2pRtxG

283◆2VB8wsVUoo 2017/11/25(土) 00:24:22.25ID:kb2pRtxG

284◆2VB8wsVUoo 2017/11/25(土) 00:24:46.55ID:kb2pRtxG

285◆2VB8wsVUoo 2017/11/25(土) 00:25:07.98ID:kb2pRtxG

286◆2VB8wsVUoo 2017/11/25(土) 00:25:29.17ID:kb2pRtxG

287◆2VB8wsVUoo 2017/11/25(土) 00:25:50.80ID:kb2pRtxG

288◆2VB8wsVUoo 2017/11/25(土) 00:26:12.77ID:kb2pRtxG

289132人目の素数さん2017/11/25(土) 01:57:31.82ID:sNF3t9bW
>>278
>argsは見れる関数と見れない関数があるのであんまり役に立たない
だったら、例が悪いよ。見られない関数の例をあげないと

> ヘルプで見ようと思ったら書いてなかった
こっちだってわざわざ下請け関数のオプションまでは書かないよ。
ヘルプを見て未解決ならソースを見てくれ。

290◆2VB8wsVUoo 2017/11/25(土) 06:17:05.29ID:kb2pRtxG

291◆2VB8wsVUoo 2017/11/25(土) 06:17:23.25ID:kb2pRtxG

292◆2VB8wsVUoo 2017/11/25(土) 06:17:39.26ID:kb2pRtxG

293◆2VB8wsVUoo 2017/11/25(土) 06:17:55.24ID:kb2pRtxG

294◆2VB8wsVUoo 2017/11/25(土) 06:18:11.86ID:kb2pRtxG

295◆2VB8wsVUoo 2017/11/25(土) 06:18:28.49ID:kb2pRtxG

296◆2VB8wsVUoo 2017/11/25(土) 06:18:44.12ID:kb2pRtxG

297◆2VB8wsVUoo 2017/11/25(土) 06:19:01.27ID:kb2pRtxG

298◆2VB8wsVUoo 2017/11/25(土) 06:19:19.56ID:kb2pRtxG

299◆2VB8wsVUoo 2017/11/25(土) 06:19:36.14ID:kb2pRtxG

300132人目の素数さん2017/11/26(日) 17:41:27.37ID:ahD9OIUm
rm(list=ls())

で R のメモリーをクリア

301132人目の素数さん2017/11/28(火) 13:00:14.06ID:kbsEO2Dt
先週長々と質問させていただいたものです。その後、試行錯誤して行った結果十解決しました。

以前は「summary」としか入れていなかったのですが、「dd$summary」と入力したところ
全体のデータが出力されました。それ以後「dd$summary$dir.area」等で出せるようになりました。
大変初期的な場面で躓いておりましたが、現在何とか進んでおります。
色々とお答えくださった方々、とても参考になりました。本当にありがとうございました。

302132人目の素数さん2017/12/10(日) 13:11:49.22ID:cG4z9akW
N(=100)回コインをなげてn(=5回)以上続けて表がでる確率。

seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
mean(replicate(10^5,seqn()))

> mean(replicate(10^5,seqn()))
[1] 0.81085

案外、高い確率になった。

303132人目の素数さん2017/12/10(日) 22:28:22.98ID:WuWuRvn7
ポアソンでやるやつだっけ

304132人目の素数さん2017/12/11(月) 07:47:07.43ID:EURogHTh
pooledVariance <- function(...) {
args = list(...)
n.args=length(args)
ss2=0
df=0
for(i in 1:n.args){
ss2 = ss2 + var(args[[i]])*(length(args[[i]])-1)
df = df + (length(args[[i]])-1)
}
ss2/df
}

effectsize <- function(y1,y2){
diff=mean(y1)-mean(y2)
var=(var(x1)*(length(x1)-1)+ var(x2)*(length(x2)-1))/(length(c(y1,y2))-2)
sd=sqrt(var)
diff/sd
}

library(effsize)
cohen.d()

305132人目の素数さん2017/12/11(月) 07:50:22.91ID:EURogHTh
>>303
単なる二項分布。
コインが5回続けて表がでたら、0.5^5 <0.05なのでイカサマコインといわれちゃいそうなんだが、
100回やってみると案外、5回表が続くので確率を計算しようと思ったが、解析的にできる頭がないので
シミュレーションしてみた。

306132人目の素数さん2017/12/11(月) 07:55:42.42ID:EURogHTh
>302は

> rbinom(100,1,0.5)
[1] 0 1 0 0 1 0 0 0 1 1 1 0 1 0 1 1 0 1 1 0 1 1 0 0 0 1 0 1 1 1 1 0
[33] 0 1 0 0 0 0 0 0 1 0 1 1 0 0 0 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 0
[65] 0 0 1 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 1 1 1 1 1 0 1 1 0 0 1 1 0 0
[97] 0 0 1 1

で5回以上1が連続するときTRUEを返す関数なのだが
もっと簡単にやれないかなぁとは思っている。

rep(1,5) %in% rbinom(100,1,0.5)は 1個ずつ評価されてTRUEが5個返ってくるだけ。

文字列にしてgrepを使うとなんとかなりそうな気がしないでもないのだけど....

307132人目の素数さん2017/12/22(金) 17:21:16.51ID:U1Bq8bdb
すまん、教えてほしいだけど
分析するために初めてRをインストールしようと思って、このスレのあるように公式サイト行ったら、esetが「JS/Redirector.NAV トロイの木馬」を検知したんだが…;
どうしたらいいだ…

308132人目の素数さん2017/12/22(金) 20:09:03.17ID:ol6XbcBE
そっちじゃなくてCRAN行け

309132人目の素数さん2017/12/22(金) 20:42:31.02ID:U1Bq8bdb
ホームページダメとかわけわかんねぇ・・・とんでもねぇな
ありがとう

310132人目の素数さん2017/12/22(金) 20:50:20.52ID:uzZ0hejh
>>307
エロサイトにアクセスしてないw?

311132人目の素数さん2017/12/22(金) 21:27:01.31ID:RNoKRtjI
>>309
CRANはここ
https://cran.r-project.org


プロジェクトの方はESETの誤検知っぽいんだよな

312132人目の素数さん2017/12/23(土) 09:28:44.66ID:n0SBd+bp
>>310
してねぇよww
totalvirusで調べると一件引っかかるし、Redirector検知だから、しょうがないね

313132人目の素数さん2017/12/24(日) 20:42:51.55ID:CT/NKMd7
col=rgb(runif(1),runif(1),runif(1),runif(1))で色指定すると

走らせるたびに色がちがっておもしろい。

hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))

314132人目の素数さん2017/12/31(日) 17:23:19.35ID:14tdpK/Y
stanやJAGSのコードでgamma関数を使おうとして
y = gamma(x)
と、書いたらエラーになった。

stanだと y=tgamma(x)、JAGSだとy=exp(loggam(x))で動作した。

315132人目の素数さん2018/01/02(火) 08:46:50.95ID:qdmBZ37O
ある大学の入学者男女の比率は1であるという帰無仮説を検定する課題が花子と太郎に課された。

花子は50人を調査できたら終了として入学者を50人をみつけて18人が女子であるという結果を得た。
帰無仮説のもとで
50人中18人が女子である確率は 0.01603475
これ以下になるのは50人中0〜18人と32〜50人が女子の場合なので
両側検定して
> sum(dbinom(c(0:18,32:50),50,0.5))
[1] 0.06490865
> binom.test(18,50,0.5)$p.value
[1] 0.06490865
で帰無仮説は棄却できないと結論した。
http://i.imgur.com/XDIp9rM.png

一方、本番と十八番が好きな太郎は一人ずつ調べて18人めの女子がみつかったところで調査を終えることにした。
18人めがみつかったのは花子と同じく50人めであった。
帰無仮説のもとで
18人がみつかるのが50人めである確率は0.005772512
これ以下になるのは23人以下50人以上番めで女子18人めがみつかった場合なので
両側検定して
pnb=dnbinom(0:999,18,0.5)
> 1 - sum(pnb[-which(pnb<=dnbinom(50-18,18,0.5))]) # < 0.05
[1] 0.02750309
http://i.imgur.com/K3T7utr.png
で帰無仮説は棄却される。

どちらの検定が正しいか、どちらも正しくないか?
検定する意図によってp値が変わるのは頻度主義統計の欠陥といえるか?

花子の横軸は裏口入学者数、太郎の横軸はサンプル数なので
サンプルでの裏口入学率を横軸にして95%信頼区間を示す。
花子の検定での信頼区間は0.36〜0.72で18/50を含む、p=0.06491
http://i.imgur.com/SeTLk8K.jpg
太郎の検定での信頼区間は0.375〜0.72で18/50を含まない、p= 0.0275
http://i.imgur.com/tNzlfxe.jpg
主観である、検定の中止の基準の差でp値や信頼区間が変化するのは変だという批判である。

316132人目の素数さん2018/01/03(水) 10:35:48.54ID:YJfyxrv+
(訂正)
ある大学の入学者男女の比率は1であるという帰無仮説を検定する課題が花子と太郎に課された。

花子は50人を調査できたら終了として入学者を50人をみつけて18人が女子であるという結果を得た。
帰無仮説のもとで
50人中18人が女子である確率は 0.01603475
これ以下になるのは50人中0〜18人と32〜50人が女子の場合なので
両側検定して
> sum(dbinom(c(0:18,32:50),50,0.5))
[1] 0.06490865
> binom.test(18,50,0.5)$p.value
[1] 0.06490865
で帰無仮説は棄却できないと結論した。
http://i.imgur.com/XDIp9rM.png

一方、十八という数字が好きな太郎は一人ずつ調べて18人めの女子がみつかったところで調査を終えることにした。
18人めがみつかったのは花子と同じく50人めであった。
帰無仮説のもとで
18人がみつかるのが50人めである確率は0.005772512
これ以下になるのは23人以下50人以上番めで女子18人めがみつかった場合なので
両側検定して
pnb=dnbinom(0:999,18,0.5)
> 1 - sum(pnb[-which(pnb<=dnbinom(50-18,18,0.5))]) # < 0.05
[1] 0.02750309
http://i.imgur.com/K3T7utr.png
で帰無仮説は棄却される。

どちらの検定が正しいか、どちらも正しくないか?
検定する意図によってp値が変わるのは頻度主義統計の欠陥といえるか?

花子の横軸は女子数、太郎の横軸はサンプル数なので
サンプルでの女子の割合を横軸にして95%信頼区間を示す。
花子の検定での信頼区間は0.36〜0.72で18/50を含む、p=0.06491
http://i.imgur.com/SeTLk8K.jpg
太郎の検定での信頼区間は0.375〜0.72で18/50を含まない、p= 0.0275
http://i.imgur.com/tNzlfxe.jpg
主観である、検定の中止の基準の差でp値や信頼区間が変化するのは変だという批判である。

317132人目の素数さん2018/01/15(月) 16:30:41.76ID:wJofbCL/
kainokousiki<-function(a,b,c){return (-b+sqrt(b^2-4*a*c))/(2*a)} #解の公式
kainokousiki(1,-5,6)

でrunすると3じゃなくて6を返すんだけど、どこが間違ってる?

318132人目の素数さん2018/01/15(月) 16:35:00.20ID:wJofbCL/
自己解決
かっこが足りなかった

319132人目の素数さん2018/01/19(金) 12:12:23.00ID:k3h1PrrK
Pythonのスレはないのか

320◆2VB8wsVUoo 2018/01/21(日) 20:11:34.35ID:oUqQkvBY

321◆2VB8wsVUoo 2018/01/21(日) 20:11:56.34ID:oUqQkvBY

322◆2VB8wsVUoo 2018/01/21(日) 20:12:22.56ID:oUqQkvBY

323◆2VB8wsVUoo 2018/01/21(日) 20:12:47.97ID:oUqQkvBY

324◆2VB8wsVUoo 2018/01/21(日) 20:13:05.14ID:oUqQkvBY

325◆2VB8wsVUoo 2018/01/21(日) 20:13:34.57ID:oUqQkvBY

326◆2VB8wsVUoo 2018/01/21(日) 20:14:02.66ID:oUqQkvBY

327◆2VB8wsVUoo 2018/01/21(日) 20:14:16.99ID:oUqQkvBY

328◆2VB8wsVUoo 2018/01/21(日) 20:14:35.47ID:oUqQkvBY

329◆2VB8wsVUoo 2018/01/21(日) 20:15:07.01ID:oUqQkvBY

330132人目の素数さん2018/02/05(月) 19:08:34.64ID:AIJT+apj
機械学習をきっかけにPythonに逆転された感じだね

331132人目の素数さん2018/02/06(火) 18:35:39.90ID:tAZA/Fp/
『Rを使った〜』だとPythonじゃないのかよって思うよね

332132人目の素数さん2018/02/06(火) 21:37:32.10ID:tUqX17n9
ずっとRだけでPython触ったこと無いけど、覚え直す価値ある?
環境構築からもう面倒なイメージ

333132人目の素数さん2018/02/06(火) 22:43:23.82ID:qQMNyZjW
Python自体は:と直後のインデントさえ気を付ければかなり簡単
3系は数が全て小数扱いなので楽

Anacondaというパッケージでインストールすれば、今流行りのJupyter Notebookという開発環境で対話的にコーディングできる(Rも使える)

334132人目の素数さん2018/02/07(水) 20:40:49.48ID:mgaw9oVw
アナコンダてのがRぽくできるのね、ありがとう
dplyrやggplot2みたいに素人でも簡単便利だといいんだけど
Pythonはオブジェクト志向ぽいしすぐ諦めそう

335132人目の素数さん2018/02/12(月) 16:19:11.21ID:NSJ4iUa4
ブラウザ環境なくなっちゃったの?

336132人目の素数さん2018/02/12(月) 17:06:56.33ID:hMO/kWPg
>>335
誤爆?
何をブラウズする環境が無くなったの?

337132人目の素数さん2018/02/13(火) 18:21:54.45ID:1hic99Cx
ブラウザでプログラミングする環境

338132人目の素数さん2018/02/13(火) 21:09:53.07ID:EFzybxrC
>>337
RStudio Serverの話?

339132人目の素数さん2018/02/13(火) 22:27:26.57ID:Qyl0hNg7
RStudio Severあるで。Dockerで使うのがいいんじゃないかな?!

340132人目の素数さん2018/02/25(日) 18:12:28.68ID:aD34K55o
統計学とウェブ解析を交えて実践的な勉強と練習を
したいのですが、おすすめな書籍やサイトはありますか。
実際に解析ツールや分析ツールを用いて
自分で分析解析してから
解答を見て適切な手順や方法、考察を
解説してくれるものが良いです。


統計学は統計検定2級の知識はありますが
ウェブ解析はテキスト読んだだけです。

341132人目の素数さん2018/02/26(月) 12:43:55.13ID:U8DvzUE2
ウェブ解析とは具体的になんですか

342132人目の素数さん2018/02/27(火) 00:54:41.17ID:hS0OJ3qQ
R studioは日本語コメント書く度にIMEが無効になったりカーソルがずれたり黒文字の予測変換が黒背景と重なって見えなくなったりと散々だわ

343132人目の素数さん2018/02/27(火) 01:06:58.09ID:GkJoHhTZ
Windows環境で使うから

344132人目の素数さん2018/02/27(火) 01:16:41.39ID:hS0OJ3qQ
会社規定なんでしかたない
UNIX環境使えるのが羨ましい

345132人目の素数さん2018/02/27(火) 01:25:59.40ID:34tpClZB
Docker for WindowsでRStudioサーバー動かせば?

346132人目の素数さん2018/02/27(火) 12:19:33.31ID:ZcuB0rsn
なぜPythonスレッドがない?

347132人目の素数さん2018/02/27(火) 14:05:09.17ID:O+8uJ5V+
板名読める?

348132人目の素数さん2018/03/04(日) 20:14:10.91ID:nWbPE6nD
>>342
IME無効になるのは俺だけじゃなかったと知ってほっとした。

349132人目の素数さん2018/03/04(日) 21:28:34.58ID:R7ZPBSuG
ストアアプリも同じ症状でるからRStudio固有の問題でなくWindows環境の不治の病だと思ってる

350132人目の素数さん2018/03/26(月) 20:57:59.83ID:WAyucZm1
Run any R code you like. There are over three thousand R packages preloaded.

https://rdrr.io/snippets/

351132人目の素数さん2018/03/26(月) 21:11:25.16ID:WAyucZm1
>>350
日本語あると動作しないが、
グラフを描いてくれるのはうれしい。


http://imagizer.imageshack.com/img923/2879/m59LVv.png

352132人目の素数さん2018/04/04(水) 07:52:09.98ID:PZp+DZN4
Fisher test検定時に
p<2.2e-16
と表示されるんですが、これより小さい値の指数桁数を正確に表記する方法教えて下さい。
例えば5.8e-35となるようにです。

353132人目の素数さん2018/04/04(水) 08:05:49.72ID:lBVECsOD
返値の中にあるp値を参照しなされ

354132人目の素数さん2018/04/04(水) 10:02:52.31ID:CRvlhZKw
fisher.test関数の返り値はlist型で、その中にp.valueという名前でp値が格納されているから$演算子を使って直接参照するか、broom::tidy関数に返り値を渡してdata.frame形式で出力してやれば見れる

355132人目の素数さん2018/04/04(水) 12:43:50.08ID:LbKgW3kd
>>352
>>354が言う直接的な参照
> fisher.test(matrix(c(1,120,130,2),2))$p.value
[1] 1.691912e-69

356132人目の素数さん2018/04/04(水) 13:06:17.36ID:PZp+DZN4
352です。よく分かりました!
ありがとうございます!

357132人目の素数さん2018/04/04(水) 14:44:33.05ID:LbKgW3kd
>>356
技術的な助言をしたけど、学術的に言えば、
p < 0.01 は全て p < 0.01 として、具体的なp値を考える意味はないと思うよ。
一部の例外的な研究分野を除いて(e.g., 遺伝統計学)。

358132人目の素数さん2018/04/04(水) 15:18:50.72ID:PZp+DZN4
はい、まさにその例外的な分野で使おうとしてます。ありがとうございます。

359◆2VB8wsVUoo 2018/04/07(土) 06:52:38.98ID:yx+HETs3

360◆2VB8wsVUoo 2018/04/07(土) 06:53:00.80ID:yx+HETs3

361◆2VB8wsVUoo 2018/04/07(土) 06:53:20.77ID:yx+HETs3

362◆2VB8wsVUoo 2018/04/07(土) 06:53:40.58ID:yx+HETs3

363◆2VB8wsVUoo 2018/04/07(土) 06:54:01.92ID:yx+HETs3

364◆2VB8wsVUoo 2018/04/07(土) 06:54:24.32ID:yx+HETs3

365◆2VB8wsVUoo 2018/04/07(土) 06:54:47.76ID:yx+HETs3

366◆2VB8wsVUoo 2018/04/07(土) 06:55:08.41ID:yx+HETs3

367◆2VB8wsVUoo 2018/04/07(土) 06:55:30.46ID:yx+HETs3

368◆2VB8wsVUoo 2018/04/07(土) 06:55:54.53ID:yx+HETs3

369132人目の素数さん2018/04/07(土) 10:24:54.55ID:maRES4NF
>>357
多重比較だと意味あるかも

370132人目の素数さん2018/04/19(木) 21:35:35.54ID:GVMUXyX9
Rのガンマ関数はいくつでオーバーフローするかやってみた。

> i=1
> while(gamma(i)!=Inf){
+ i=i+1
+ }
Warning message:
In gamma(i) : value out of range in 'gammafn'
> i
[1] 172
> gamma(172)
[1] Inf
Warning message:
value out of range in 'gammafn'
> gamma(171)
[1] 7.257416e+306

371132人目の素数さん2018/04/29(日) 00:21:51.67ID:5dW+xNwa
matplot()で折れ線グラフ描いたときに、X軸をカテゴリで示したいのですが、
可能でしょうか?

例えばtemp <- c("0時間","8時間","24時間","48時間")として、
matplot()の引数にtempをとるやり方です。
他にもやり方あれば教えてください。

372132人目の素数さん2018/04/29(日) 01:55:59.75ID:PngyHOQZ
>>371
matplot(..., xaxt="n")
axis(1, at=seq(along=temp), lab=temp)

373132人目の素数さん2018/04/30(月) 16:50:51.65ID:t3vhzyao
>>372
遅くなりましたがありがとうございました。
できました!

374◆2VB8wsVUoo 2018/04/30(月) 23:55:02.74ID:y1TqbRSE

375◆2VB8wsVUoo 2018/04/30(月) 23:55:22.50ID:y1TqbRSE

376◆2VB8wsVUoo 2018/04/30(月) 23:55:43.19ID:y1TqbRSE

377◆2VB8wsVUoo 2018/04/30(月) 23:56:02.89ID:y1TqbRSE

378◆2VB8wsVUoo 2018/04/30(月) 23:56:23.17ID:y1TqbRSE

379◆2VB8wsVUoo 2018/04/30(月) 23:56:40.16ID:y1TqbRSE

380◆2VB8wsVUoo 2018/04/30(月) 23:56:53.49ID:y1TqbRSE

381◆2VB8wsVUoo 2018/04/30(月) 23:57:13.43ID:y1TqbRSE

382◆2VB8wsVUoo 2018/04/30(月) 23:57:32.05ID:y1TqbRSE

383◆2VB8wsVUoo 2018/04/30(月) 23:57:53.02ID:y1TqbRSE

384132人目の素数さん2018/05/01(火) 18:57:32.34ID:iUBwAKWd
特定の長方形の中に複数の長方形を最小面積で敷き詰める平面充填に関するパッケージってありませんかね

385132人目の素数さん2018/05/01(火) 19:35:44.26ID:iUBwAKWd
やっぱ自分でアルゴリズムを書いてみます

386132人目の素数さん2018/05/02(水) 00:09:39.04ID:xpXz1bJY
Package`deldir'
これどうなの?

387132人目の素数さん2018/05/06(日) 22:15:43.42ID:BK1CxH7U
# jonckheereテストを書いてみた

jonckheere <- function(L,
alternative = c("two.sided", "increasing", "decreasing"),
cat=TRUE){
# L : list of vectors A1,A2,...,Ak, with assumed tendency
How.Many.Greater.Pairs <- function(A,B){ # How many pairs of A[i] > B[j], count as 0.5 when equal,
n.a = length(A)
n.b = length(B)
how.many.greater.pairs = 0
for(i in 1:n.a){
for(j in 1:n.b){
how.many.greater.pairs = how.many.greater.pairs+ifelse(A[i]==B[j],0.5,A[i]>B[j])
}
}
return(how.many.greater.pairs)
}
Sum.of.Greater.Pairs <- function(L){ #L=list(A1,,,,Ak),A1 < A2 < A3,..,< Ak : vector
k = length(L)
comb = combn(1:k,2) # possible combinaition of pairs to compare
n.comb = ncol(comb) # how many combinations
J = 0 # sum of greater pairs
for(i in 1:n.comb){
J = J + How.Many.Greater.Pairs(L[[comb[1,i]]],L[[comb[2,i]]])
}
return(J)
}
J = Sum.of.Greater.Pairs(L)
n = sapply(L,length)
N = sum(n)
EJ = (N^2-sum(n^2))/4
VJ = (N^2*(2*N+3)-sum(n^2*(2*n+3)))/72
Z = (J-EJ)/sqrt(VJ)
alternative = match.arg(alternative)
p.value = switch(alternative, 'two.sided' = 2 * min(pnorm(Z), pnorm(-Z), 0.5),
'increasing' = pnorm(Z),
'decreasing' = pnorm(-Z))
if(cat){
cat( 'p.value = ', p.value,'\n')
cat('alternative hypothesis: ' ,alternative,'\n')
}
invisible(p.value)
}

388132人目の素数さん2018/05/22(火) 20:34:14.53ID:8iM7y57s
代入は「<-」派?「=」派?

389132人目の素数さん2018/05/22(火) 21:16:34.39ID:iB1pjrmI
>>388
実態調査か何か?
<-と=は挙動が違う場合があるので、使い分けていますが、
代入はどっちかと問われたら、無論 <- または ->

なお、
> 1 -> x
これはエラーにならないけど、
> 1 = x
1 = x でエラー: 代入の左辺が不正 (do_set) です
これはエラー

390132人目の素数さん2018/05/23(水) 11:17:56.08ID:OSJ/4EBd
>>389
俺は基本=派。

関数の定義は
z.test <- function(x,n=16,sigma=1){
z=sqrt(n)*mean(x)/sigma
2*pnorm(abs(z),lower=FALSE)
}
と書いている。

391132人目の素数さん2018/05/23(水) 11:24:57.67ID:MGQGuwX9
>1 = x でエラー

当たり前
そんな使い方なんてするかよ

他言語と同じく=一文字の方がすっきりしてイイ

392132人目の素数さん2018/05/24(木) 14:28:38.60ID:ExPgBsbL
こういうのが紛らわしいから、俺は = 推奨。

x <- 1
if(x <- 1) print('YES')
if(x < -1) print('YES')

393132人目の素数さん2018/05/24(木) 14:36:00.11ID:EQ5K0CF7
だよねぇ
<-良くない

394132人目の素数さん2018/05/25(金) 08:05:58.04ID:ZHt2t+40
やったことなかったので関数の初期値設定に<-を使うとどうなるかやってみた。
まず、= の場合
> z.test <- function(x,n=16,sigma=1){
+ z=sqrt(n)*mean(x)/sigma
+ 2*pnorm(abs(z),lower=FALSE)
+ }
> z.test(1:3)
[1] 1.244192e-15

<- で初期値設定すると、エラー

> z0.test <- function(x,n<-16,sigma<-1){
Error: unexpected assignment in "z0.test <- function(x,n<-"
> z=sqrt(n)*mean(x)/sigma
Error in mean(x) : object 'x' not found
> 2*pnorm(abs(z),lower=FALSE)
Error in pnorm(abs(z), lower = FALSE) : object 'z' not found
> }

395132人目の素数さん2018/05/25(金) 08:11:48.02ID:ZHt2t+40
俺は見栄えがいいと思って関数定義には<-を使っているけど
= でも通常に動作する。
> z.test = function(x,n=16,sigma=1){
+ z=sqrt(n)*mean(x)/sigma
+ 2*pnorm(abs(z),lower=FALSE)
+ }
> z.test(1:3)
[1] 1.244192e-15

<- 推奨の人に聞きたいのだけど
<- でないと動作しないってことあるのだろうか?

396132人目の素数さん2018/05/25(金) 09:41:36.92ID:GCdXSt8a
<-は、打つのが
=より面倒

うっとおしい

397132人目の素数さん2018/05/25(金) 20:19:41.68ID:dXqzteX1
RStudioつかってりゃ[Alt+-]で簡単入力

それと引数の指定は代入じゃないと思うんだが感覚が違うよかな?

398132人目の素数さん2018/05/25(金) 21:33:31.44ID:QNMt6z2O
>>396
essだとアンダースコアを入れるを勝手に「<-」になる。
>>394
それは代入ではなく、関数の規定値の設定。
規定値の設定は「=」と決められているので、「<-」はアウト。
ただし、関数を実行するときには「<-」を使うことができる。
> mean(x<-1:10)
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10

399132人目の素数さん2018/05/25(金) 22:26:21.26ID:Cz/+g7h9
=に拘り続けて、<-に変える瞬間の君たちがたまらないよ

400132人目の素数さん2018/05/28(月) 16:26:17.30ID:d/09kgU6
>>398
ありがとうございます。
そんな使い方ができたのですね。

こんなのしか知りませんでした。

> f <- function() {
+ x<<-1:10
+ mean (x)
+ }
> f()
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10

401132人目の素数さん2018/05/28(月) 20:56:51.59ID:Osxttqv4
一画面しかグラフがないのに

Hit <Return> to see next plot:

と出るの鬱陶しいな、と思っていたら

par(ask=FALSE)

と設定しておけばいいんだな。

402132人目の素数さん2018/05/29(火) 18:08:52.29ID:ZCWRCH5Y
同名のdataがあるのでパッケージを::で指定してもうまくいかなかった。

data(netmeta::parkinson)

でなくて

data(parkinson, package = 'netmeta')

とするんだな。

403132人目の素数さん2018/05/30(水) 06:16:03.29ID:fLd3NENr
hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
だと、1色だけど

hist(rnorm(100),col=sample(colours(),(sample(1:10,1))))
なら1〜10色で表示される。

404132人目の素数さん2018/05/31(木) 21:37:03.79ID:BfVjjX7C
runif(1)なんて・・・

405132人目の素数さん2018/05/31(木) 23:16:15.39ID:grs1zCKo
> hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
何色に増やしてもいいがな。
何をやりたいのかな?
hist(rnorm(5000),col=apply(matrix(runif(80),4), 2, function(x){rgb(x[1],x[2],x[3],x[4])}))

406132人目の素数さん2018/06/01(金) 21:54:20.67ID:Ef899k/0
:::てどういう時に使うんだろ?
ソース読みたくて
library(BayesFactor)
ttestBF_indepSample
では表示されなかったが、
library(BayesFactor)
BayesFactor:::ttestBF_indepSample
だと出てきた。

407132人目の素数さん2018/06/04(月) 13:04:13.10ID:3QOAvZa7
Stanも覚えなきゃいけないのかよーくそがー

408132人目の素数さん2018/06/04(月) 13:05:22.57ID:3QOAvZa7
::なら名前空間カブールよけよけの術だよね

409132人目の素数さん2018/06/04(月) 14:24:00.43ID:d7l4LqX9

410132人目の素数さん2018/06/05(火) 12:44:14.88ID:LfwL/Kok
>>408
:::と3個なんだよなぁ

411132人目の素数さん2018/06/05(火) 13:06:03.27ID:CHJNgrVq
>>410
グダグダ言っていないで、さっさとヘルプを参照すればいいじゃないか?
?':::'

412132人目の素数さん2018/06/05(火) 19:45:28.11ID:LfwL/Kok
>>411
ありがとう。

stats:::t.test.default

でt検定のソースが見られた

413132人目の素数さん2018/06/11(月) 00:14:18.19ID:K1zGYlRd
関数を値で上書きしてしまった場合に元に戻す方法を教えて頂きたく。

例)
rm <- 1
としてしまった場合にrm関数を元に戻したいです。


よろしくお願いします。

414132人目の素数さん2018/06/11(月) 07:21:35.69ID:+ayC4RGv
>>413
直後ならRstudioでctrl+zで元に戻す。
再起動してたらこの方法は無理だろな

415132人目の素数さん2018/06/11(月) 08:24:34.24ID:2ZHCUqeW
>>413
自分で作ったのではない関数の rm だったら rm(rm) で数値ベクトルオブジェクトが消えて、関数の rm が見えるようになる。

416132人目の素数さん2018/06/11(月) 15:01:37.34ID:mlFyU0v4
>>413
> rm <- 1
> ls()
[1] "rm"
> rm(rm)
> ls()
character(0)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]

深く考えなくてもrm(rm) で元に戻るけど。
rmが数値ではなくて自作関数だったら、少々ややこしいけど、
base::rm()で大丈夫だろう

417132人目の素数さん2018/06/11(月) 17:04:35.54ID:Z+okZT62
>>413
既に説明されてるけど上書きじゃなくて別オブジェクト扱いになるからrmで消すか明示的にパッケージを::演算子で関数を指定してやればおけ

自作関数だと上書きされちゃうけど

418132人目の素数さん2018/06/11(月) 19:33:14.14ID:rfEFkwvW
rmが自作関数だと戻せないのでは?

> rm = function(x) asin(sqrt(x))
> rm(1)
[1] 1.570796
> rm
function(x) asin(sqrt(x))
> rm = 1
> rm
[1] 1
> rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
inherits = FALSE)

419132人目の素数さん2018/06/11(月) 19:51:41.60ID:Z+okZT62
同じ名前空間にあるなら上書きされちゃうよ

.Gloalenv(だったかな?)に自作関数があり、そこに上書きして変数にしてるから消したら消えるだけで戻せない

420132人目の素数さん2018/06/11(月) 19:55:10.50ID:Z+okZT62
ちょっと古いけど、この辺りを読むとなぜそうなるかが分かると思う

(Rの)環境問題について その1。
ttps://qiita.com/kohske/items/325bdf48f4f4885a86f1

421132人目の素数さん2018/06/11(月) 20:04:41.37ID:mlFyU0v4
>>418
戻せます。
> rm <- function(x) asin(sqrt(x))
> rm
function(x) asin(sqrt(x))
> base::rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]

422132人目の素数さん2018/06/11(月) 20:06:37.76ID:mlFyU0v4
>>419
>.Gloalenv
case sensitiveなRでそれはないのでは?
> help(".GlobalEnv")

423132人目の素数さん2018/06/11(月) 20:19:14.99ID:Z+okZT62
フォローありがとう

4244132018/06/11(月) 21:00:49.78ID:K1zGYlRd
>>414-422
みなさまご返信ありがとうございました。
例示したものがbase::rmだったので自作関数のことまでは想像だにしませんでした…(がとても良い勉強になりました。)

自環境でも上書きした関数を元に戻すことが出来ることを確認しました。

425132人目の素数さん2018/06/12(火) 07:42:40.88ID:9RfEtlLW
>>421
戻したいのが自作関数
rm = function(x) asin(sqrt(x))
だったら上書きされて無理じゃないの?
ゴミ箱から消去したファイルを復活させるソフトもあるから
PCに詳しければ可能かもしれないけれど
Rの機能としては復活させるのは無理と思う。

426132人目の素数さん2018/06/12(火) 11:00:45.94ID:cD/mQB+6
自作のrm関数はbaseパッケージとは別に定義されるから
base::rm(rm)で自作の方のrmは消えてbase::rmが「復活」します。

427132人目の素数さん2018/06/12(火) 12:07:52.53ID:cD/mQB+6
ああすまん
自作の方を上書きした場合は、そっちの復活は無理ですね。

428132人目の素数さん2018/06/12(火) 14:36:07.63ID:VF7OcVBc
>>425
ちょっと意味がわからない。
自作の関数なら定義を再実行するだけなのに、
消える消えないって何が問題なんだ?

429132人目の素数さん2018/06/13(水) 14:39:27.58ID:ZXkz+qgi
>>424
関数だけでなく既定値でも同じ。

> pi
[1] 3.141593
> pi=2
> pi
[1] 2
> rm(pi)
> pi
[1] 3.14159

430132人目の素数さん2018/06/14(木) 09:27:55.98ID:mxBGyFKT
共同ツール 1
https://seleck.cc/685

https://trello.com/
ボードのメニュー → Power-Upsから拡張可能 Slack DropBoxなど
Trello Chrome拡張機能 elegant
ttp://www.kikakulabo.com/service-eft/
trelloのオープンソースあり

共同ツール 2
https://www.google.com/intl/ja_jp/sheets/about/

共同ツール 3
https://slack.com/intl/ja-jp
https://www.dropbox.com/ja/
https://bitbucket.org/
https://ja.atlassian.com/software/sourcetree
https://www.sketchapp.com/
ttp://photoshopvip.net/103903
ttps://goodpatch.com/blog/sketch-plugins/

Trello Chrome拡張機能プラグイン集
https://chrome.google.com/webstore/search/trello?_category=extensions

Slackプラグイン集
https://slack.com/apps

Sketchプラグイン集
https://sketchapp.com/extensions/plugins/
https://supernova.studio/

431132人目の素数さん2018/06/16(土) 19:29:24.89ID:10mr4q2e
はっとデシャング

432132人目の素数さん2018/06/17(日) 19:06:33.58ID:OYjqtCQI
>>306
知恵袋に漸化式があったので計算スクリプトを書いてみた。
表がでる確率pのコインをN回投げてK回以上表が続く確率。

# N: flips
# K: least sequential head
# p: probability of head
seqNp <- function(N=100,K=5,p=0.5){
if(N==K) return(p^K)
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q for any i

for(i in K:(N-1)){ # recursive formula
a[i+1]=0 # a[i+1]=q/p*(a[i]+a[i-1]+a[i-2]+...+a[i-(K-1)])
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}

P0=numeric(N) # P0[n] : probability of ending with 0 head when flipped n times
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n

MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0
# MP[k,n] = Pk[n] : probability of ending with k head when flipped n times
for(k in 2:K){ # Pk(n+1)=p*P(k-1)(n)
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
}
ret=1-apply(MP,1,sum)

ret[N]
}

433132人目の素数さん2018/06/18(月) 09:16:00.35ID:JdABdagZ
>>306
この例で、答えとして何が返ってきて欲しいのか?
2で合ってる?

434132人目の素数さん2018/06/18(月) 21:23:58.67ID:fHaLluDH
>>433
1が5回以上出現しているからTRUE(数字なら1)が返り値。

435132人目の素数さん2018/06/19(火) 10:20:04.23ID:fcA67BpN
>>434
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.

chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}

436132人目の素数さん2018/06/19(火) 21:57:17.85ID:QuKrVrE8
>>435
whichとdiffを用いての解法のお返事ありがとうございます。

実行時間を>302のスクリプトと比べてみました。

run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.

chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}


# N(=100)回コインをなげてn(=5回)以上続けて表がでる確率。
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}

> system.time(mean(replicate(10^5,run.check(rbinom(100,1,0.5)))))
user system elapsed
17.56 0.04 17.68
> system.time(mean(replicate(10^5,seqn())))
user system elapsed
9.74 0.07 9.78

後者の方が速いのは1がn個連続したら、TRUEを返して終了するので以後のチェックはしないためだろうと思います。

diff や anyの使い方が勉強になりました。

時間を割いてスクリプトを作成していただいてありがとうございます。<(_ _)>

437132人目の素数さん2018/06/21(木) 10:12:41.96ID:Ze2kEGUX
Rstudioで日本語の入力ができないのですが、解決する方法知ってる方いませんか?
具体的には、IME自体をオンにすることができません。
日本語の表示やコピペは問題なくできます。

環境
Linux - Utubntu
RStudio Version 1.0.143

438132人目の素数さん2018/06/21(木) 10:26:53.84ID:Fds/4hNR
RStudioのQtが古いのでパッチ当てないと日本語入力できないのが現状ですわ
パッチはUbuntu 16.04.LTS用にしかないので他のバージョンならRStudioサーバー使った方が早いかと

なお、パッチはRStudio 日本語とかでグクってみて

439132人目の素数さん2018/06/23(土) 17:38:43.39ID:aZCZP6wm
windows版のRStudioでもたまに日本語入力ができなくなる。
別のソフトでIMEのON/OFFをやってから戻ると直るけど。

440132人目の素数さん2018/06/23(土) 20:18:55.21ID:4XxgOFIA
>>439
うちも同じ。

441132人目の素数さん2018/06/23(土) 21:30:47.03ID:JcH4CpHc
Windowsはストアアプリでも同じ症状でるから持病なんじゃないかと思う
Ubuntuはfcitxのパッチ当てとけばWindowsのような症状は出ないので快適

442132人目の素数さん2018/06/25(月) 13:23:09.19ID:rmZ6BvQE
これをやってみた。
http://img-cdn.jg.jugem.jp/d75/1405116/20140609_1105745.jpg

スクリプトは
http://imagizer.imageshack.com/img921/8020/t1VcqA.jpg
小数点以下500桁でみると19個素数があった。

> vip=Vectorize(is.prime)
> (idx=which (vip(N)))
[1] 100 124 150 172 183 202 215 219 255 296 310 314 323 346 400 417 439 441 474
> cat(substring (e,idx[1],idx[1]+9))
7427466391

443132人目の素数さん2018/06/26(火) 07:08:45.35ID:Na/Ih9Bj
問題

99人の囚人がいます。彼らの頭に1〜100までのナンバーカードが貼りつけられた帽子をランダムにかぶせます。
他人の帽子は見ることができても、自分の帽子は見ることができません。
帽子の数は全部で100なので、一つ使われずに余ります。
そのナンバーは囚人達にはわからないようにしておきます。
この状況で、囚人たちに一斉に自分のナンバーを宣言させて、全員が正解だったら釈放するという賭けをします。
囚人たちには帽子をかぶせられる前に相談タイムが設けられています。
どういう戦略を取れば、助かる確率を最も高くできるでしょうか?
http://000013.blogspot.com/2010/12/99.html

解答を読んでも数理が理解できないが
確率が0.5になるのはシミュレーションできた。

# http://000013.blogspot.com/2010/12/99.html
inversion <- function(x){
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret) # inversion number
}
is.even= function(x) !inversion(x)%%2 # is inverion number even?

prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx])
Y=numeric(n-1)
for (i in 1:(n-1)){ # select as even permutation
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # two numbers unseen for i-th prisoner
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2])
}
all(X==Y)
}
mean(replicate(1e3,prisoner99()))

444132人目の素数さん2018/07/11(水) 16:20:29.88ID:zFHa28EV
>>437
遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ

445132人目の素数さん2018/07/22(日) 11:22:38.35ID:Ott8rTSz
覆面算 RED + WHITE = COLOR

どうもCでやったらうまくいかないのでRでやってみた。

# RED + WHITE = COLOR
x=c('R','E','D','W','H','I','T','E','C','O','L','O','R')
unique(x)
redwhite <- function(R,E,D,W,H,I,T,C,O,L){
red=10^(2:0)
white=10^(4:0)
color=10^(4:0)
sum(red*c(R,E,D))+sum(white*c(W,H,I,T,E)) - sum(color*c(C,O,L,O,R))
}

x=unique(x) ; x
REDWHITECOLOR <- function(x){
R=x[1]
E=x[2]
D=x[3]
W=x[4]
H=x[5]
I=x[6]
T=x[7]
C=x[8]
O=x[9]
L=x[10]

cat(paste('RED = ',R,E,D,
' WHITE = ',W,H,I,T,E,
' COLOR = ',C,O,L,O,R),'\n')
}
REDWHITECOLOR(unique(x))

library(gtools)
perm=permutations(n=10,r=10,v=0:9)
perm=perm[perm[,1]!=0&perm[,4]!=0&perm[,8]!=0,] # R!=0,W!=0,C!=0
head(perm) ; tail(perm)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=redwhite(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6],perm[i,7],perm[i,8],perm[i,9],perm[i,10])
}
hist(re)
(indx=which(re==0))
REDWHITECOLOR(perm[indx,])

> REDWHITECOLOR(perm[indx,])
RED = 5 8 7 WHITE = 3 9 6 1 8 COLOR = 4 0 2 0 5

446132人目の素数さん2018/07/22(日) 12:50:02.38ID:Ott8rTSz
# AAB
# × CD
# ------
# CCE
# BFAA
# ------
# BEBDE
f <- function(A,B,C,D,E,F){
(A*100+A*10+B)*D==C*100+C*10+E &(A*100+A*10+B)*C==B*1000+F*100+A*10+A &
(C*100+C*10+E)+(B*1000+F*100+A*10+A)*10==B*10000+E*1000+B*100+D*10+E
}
library(gtools)
perm=permutations(n=10,r=6,v=0:9)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=f(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6])
}
indx=which(re==TRUE)
perm[indx,]

447132人目の素数さん2018/07/23(月) 11:09:01.32ID:SDW/E7TY
大変初歩的な質問で恐縮ですが、ご教示いただければ幸いです。

PERT分布を学んでいるのですが、ベータ分布のスケーリングの仕方が分かりません。
最小値a、最頻値b、最大値c として、a, b, c から求めたパラメータをa1、a2 とします。
このとき、
dbeta(x, a1, a2)*(c-a)+a
とすると、当然ながら横軸の値の範囲は0から1で、縦軸の値のみ大きくなってしまいます。

横軸の値の範囲がaからcで、縦軸の値は確率密度のままにするには、どうしたらよいのでしょうか?

448132人目の素数さん2018/07/23(月) 12:15:25.71ID:llIMc79p
え?aだけ右に動かしたいならこれでいいのでは?
dbeta(x-a, a1, a2)

449132人目の素数さん2018/07/23(月) 18:17:20.65ID:SDW/E7TY
>>448

447です。ご教示どうもありがとうございました。

a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
mu <- (a+4*b+c)/6 # PERT分布の平均値
a1 <- b*(mu-a)/(c-a) # パラメータ a1
a2 <- b*(c-mu)/(c-a) # パラメータ a2
pert1 <- function(x) dbeta(x,a1,a2) # ベータ分布布の計算


まではできたのですが、ご教示いただいた方法で
pert1 <- function(x) dbeta(x-a ,a1,a2)
とすると、x の範囲が0から1で、yが0の一様分布のようなグラフになってしまいます。

450132人目の素数さん2018/07/23(月) 20:05:23.59ID:5MGL8f5L
4パラべーたへの変換は下記かな?
https://en.wikipedia.org/wiki/Beta_distribution#Four_parameters_2
a1とa2のパラメータ推定って合ってるんだっけ?僕も正直わかんない

やりたいことはこれだと思う
install.packages("mc2d")してからやってね
----------------------------------------------------
rm(list=ls())
library(mc2d)

a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値

pert4<-function(x) dpert(x,min=a,mode=b,max=c,shape=4)
x1<-40:70
plot(x1,pert4(x1))

451132人目の素数さん2018/07/23(月) 21:05:56.52ID:+HoGXx9d
>>450
447です。

450様、まさにこれ、これです。どうもありがとうございました。
おかげさまで大変助かりました。

448様にも、ご助言下さいましたこt、厚く御礼申し上げます。

452132人目の素数さん2018/07/27(金) 15:50:12.19ID:E7dz8jqT
lmtestのパッケージをインストールしようとすると不正なマルチバイト文字がありますと出てくるのですがどうすればいいですか?

453132人目の素数さん2018/07/27(金) 17:33:01.36ID:ufA3ZDTu
>>452
Windows 10 (Rの日本語ヘルプはインストールしていない)+ R studioだけど、そのメッセージ出ないな。

454132人目の素数さん2018/07/27(金) 23:29:04.58ID:klPq25x/
ユーザー名が日本語かな

455132人目の素数さん2018/07/28(土) 12:56:05.38ID:lhV9Awbf
職場のSPSS信者たちにRを認めてもらうにはどうしたら良いだろう
「安かろう悪かろう」「俺が知らないものは三流」的な考え方の人が多くて肩身が狭い…

456132人目の素数さん2018/07/28(土) 17:53:14.05ID:XHIBT/Gx
>>455
>>452のようなエラーが出たら「無料のものはやっぱりダメね」
有償(SPSS) →信頼できる
ボランティア→いい加減
と言う思考らしい

457132人目の素数さん2018/07/28(土) 18:33:45.13ID:l8Bhwsl9
>>454
ユーザー名はアルファベット1文字です
studioインストールしてないんですけど、しないとダメですか?

458132人目の素数さん2018/07/28(土) 21:48:13.01ID:hNHqu0EH
他のパッケージはインストールできるの?

459132人目の素数さん2018/07/29(日) 07:46:06.57ID:u2P99O/7
linuxなら環境変数LANGをCにすれば解決する問題だからwindowsでも似たような対処すれば行けそう。

460132人目の素数さん2018/07/29(日) 19:38:42.16ID:hjP1s1BE
これかな?

> Sys.setlocale(locale="C") # locale の変更(USの設定にしないと表示が変になる)
[1] "C"
元に戻すには
> Sys.setlocale(locale="Japanese_Japan.932")

461132人目の素数さん2018/08/02(木) 20:59:18.66ID:DvtHDake
すみません、教えて下さい。

フィッシャーの正確確率検定で、
適合度(一様性)の検定ってできますか?
できるのは独立性の検定のみ?

Rでいろいろ試したみたんですが、
2群ないとエラーになる。。

462132人目の素数さん2018/08/02(木) 21:16:50.41ID:GWs8Hfcd
適合度ならカイ二乗検定でやればいいのでは?

463132人目の素数さん2018/08/02(木) 21:52:11.65ID:DvtHDake
>>462
レス感謝です。

そうなんですが、
期待度数が5未満になってしまった場合は、
カイ二乗検定は向かないんですよね?

464132人目の素数さん2018/08/02(木) 21:54:12.25ID:82PNyNf0
それで他の方法を探しています

465132人目の素数さん2018/08/02(木) 23:02:57.38ID:clch7kxN
>>463
こういうので
r1=5;r2=4;n1=10;n2=12
prop.test(c(r1,r2),c(n1,n2),correct=TRUE)
chisq.test(matrix(c(r1,r2,n1-r1,n2-r2),2,byrow=TRUE),correct=TRUE)

警告がでるという話かな?
Warning message:
In chisq.test(matrix(c(r1, r2, n1 - r1, n2 - r2), 2, byrow = TRUE), :
Chi-squared approximation may be incorrect

466132人目の素数さん2018/08/03(金) 09:04:31.58ID:VGe+pklE
>>463
適合度を見るのに関係ないだろ。
> cnt
[1] 1 0 2 4 3 4 12 6 12 11 21 22 21 37 40 30 59 44 49 47 38 48 36 43 38
[26] 46 33 34 30 21 26 33 17 13 14 12 8 16 7 7 9 6 6 5 3 8 4 4 4 1
[51] 0 2 1 0 0 0 0 2
> prb
[1] 0.0002947255 0.0006464052 0.0012538873 0.0022037831 0.0035720843
[6] 0.0054114364 0.0077413658 0.0105430123 0.0137588816 0.0172972083
[11] 0.0210398698 0.0248524558 0.0285950600 0.0321325300 0.0353432214
[16] 0.0381256526 0.0404028045 0.0421240952 0.0432652773 0.0438266419
[21] 0.0438299769 0.0433147339 0.0423338191 0.0409493611 0.0392287274
[26] 0.0372409836 0.0350539129 0.0327316483 0.0303329196 0.0279098757
[31] 0.0255074215 0.0231629878 0.0209066521 0.0187615256 0.0167443293
[36] 0.0148660914 0.0131329065 0.0115467110 0.0101060386 0.0088067297
[41] 0.0076425765 0.0066058952 0.0056880193 0.0048797146 0.0041715201
[46] 0.0035540195 0.0030180499 0.0025548579 0.0021562071 0.0018144483
[51] 0.0015225575 0.0012741483 0.0010634658 0.0008853653 0.0007352804
[56] 0.0006091856 0.0005035529 0.0004153084
> chisq.test(cnt, p=prb, rescale.p=TRUE, simulate.p.value=TRUE)

Chi-squared test for given probabilities with simulated p-value (based
on 2000 replicates)

data: cnt
X-squared = 65.504, df = NA, p-value = 0.2189
5未満の数があっても関係がない

467132人目の素数さん2018/08/03(金) 10:29:47.68ID:jRZN4rrt
>>466
警告が出るのは独立性検定の場合でしたね。

468132人目の素数さん2018/08/04(土) 02:59:15.87ID:8dNre62F
463です。
カイ二乗検定で、期待度が5未満が多いのは不適切なのかと思い込んでおりましたが、
それは独立性の検定の場合のことであり、
適合度(一様性)の検定の場合は、
気にしなくて良かったのか。。
ならカイ二乗検定でやります。

いろいろ教えて下さり、ありがとうございます。

469132人目の素数さん2018/08/04(土) 19:49:50.94ID:3IVOmTrE
なんでお前らSPSS使わないの?

http://www.cdleidyba.lt/file/spss_old_1.png

470132人目の素数さん2018/08/04(土) 20:37:20.39ID:kYh4t8ZN
>>469
Stanとかに対応してんの?

471132人目の素数さん2018/08/04(土) 22:03:56.61ID:kTWj7ufe
>>469
化石みたいなバージョンやな

472132人目の素数さん2018/08/05(日) 00:56:32.45ID:biDT5wEJ
>>469
いまSPSSを使う理由って何なの?煽りじゃなくマジで

473132人目の素数さん2018/08/05(日) 14:27:54.32ID:dSJLM/da
maple使えよ

474132人目の素数さん2018/08/05(日) 19:10:34.43ID:mULz5b3r
>>472
周囲のspssユーザを見ると「指導教員がspssしか使えないので、
自分もspssしか使えない、
新しい(学術系)ソフトは誰かに丁寧に根気よく指導してみらえないと無理」、
これの再生産。

475132人目の素数さん2018/08/05(日) 19:11:30.48ID:zaJ/WJYb
>>469

研究室に予算がないから

476132人目の素数さん2018/08/05(日) 19:12:42.57ID:mULz5b3r
>>475
予算がなかったら、Rを使え

477132人目の素数さん2018/08/05(日) 19:13:59.36ID:mULz5b3r
>>469>>475を間違えたorz

478132人目の素数さん2018/08/05(日) 19:33:14.43ID:mTZcaJ/n
予算がないならケーキを食べればいいのに

479132人目の素数さん2018/08/08(水) 13:11:20.17ID:7/3l/PxD
再帰呼び出しに関数名の代わりにRecall使うと処理に要する時間が増えることに気づいた。

フィボナッチ数列の再帰呼び出し
# 関数名での呼び出し
fibo <- function(n){
if(n==1|n==2) return(1)
else fibo(n-1)+fibo(n-2)
}

# Recallを使っての呼び出し
fiboR <- function(n){
if(n==1|n==2) return(1)
else Recall(n-1)+Recall(n-2)
}


> system.time(fibo(30))
user system elapsed
4.27 0.02 4.36
> system.time(fiboR(30))
user system elapsed
6.65 0.00 6.70

480132人目の素数さん2018/08/09(木) 17:47:22.65ID:l1M8GWWv
cumsumを再帰関数で書いてみた。

何度も試行錯誤

# cumsum with recursive call
cumsumR <- function(v){
res=numeric(length(v))
cumsumR_sub <- function(v,res,i){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
cumsumR_sub(v,res,1)
}
cumsumR(1:10)

とりえあず、動作した。
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55

481132人目の素数さん2018/08/09(木) 18:31:24.68ID:l1M8GWWv
簡略化できた

cumsumR <- function(v,res=NULL,i=1){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}

> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
>

482132人目の素数さん2018/08/10(金) 18:41:42.48ID:Hlm8Oe3x
ifelseの動作がどうも納得できないでご教示いただきたいのですが。
これってifelseの仕様でしょうか?バグでしょうか?

> x=2:1
> if(x[1]<=x[2]) x else x[2:1]
[1] 1 2
> ifelse(x[1]<=x[2],x,x[2:1])
[1] 1

483132人目の素数さん2018/08/10(金) 18:55:39.89ID:vdkgfABT
>>482
ifelse() はベクトル演算できる関数。
だから、戻り値の要素数は、test部分の要素数に一致する。
testが1つなので、x[2:1]の1つ目の要素が返される。

4844832018/08/10(金) 18:58:51.59ID:vdkgfABT
testの要素を2つにすれば、出力も2つになる。
> ifelse(c(x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2
さらに3つにすると、helpに書いている通り、yesやno部分は再利用される。
> ifelse(c(x[1] <= x[2], x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2 1

485132人目の素数さん2018/08/10(金) 23:06:17.77ID:Hlm8Oe3x
>>483
ありがとうございました。

バグではなくてそういう仕様だったのですね。

こういう使い方ができるということがわかって勉強になりました。

> x=2:1
> ifelse(c(x[1]!=x[2],x[1]==x[2]),1:2,3:4)
[1] 1 4

486132人目の素数さん2018/08/10(金) 23:17:14.10ID:VlOUWHrC
教わったので早速、

ifelse() はベクトル演算できると、再利用されるの動作確認。

> ifelse(c(TRUE,TRUE,FALSE,FALSE,TRUE,FALSE,FALSE),1:2,11:15)
[1] 1 2 13 14 1 11 12

487132人目の素数さん2018/08/10(金) 23:45:24.37ID:VlOUWHrC
2種類の文字の変換もifelseでできるんだなぁ。
> mat
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 0
[3,] 1 1 0
[4,] 1 1 0
[5,] 1 1 0
[6,] 1 1 0
> print(apply(mat,2, function(x) ifelse(x==1,'真','偽')),quote=F)
[,1] [,2] [,3]
[1,] 真 真 真
[2,] 真 真 偽
[3,] 真 真 偽
[4,] 真 真 偽
[5,] 真 真 偽
[6,] 真 真 偽

488132人目の素数さん2018/08/11(土) 00:28:09.90ID:OesSEWnz
行列も次元付きのベクトルだからもっと簡略化できた。

> (mat=matrix(sample(0:1,12,rep=T),3,4))
[,1] [,2] [,3] [,4]
[1,] 1 1 0 1
[2,] 0 1 0 0
[3,] 1 0 0 1
> print(ifelse(mat==1,'Right','Wrong'),quote=F)
[,1] [,2] [,3] [,4]
[1,] Right Right Wrong Right
[2,] Wrong Right Wrong Wrong
[3,] Right Wrong Wrong Right
>

489132人目の素数さん2018/08/11(土) 18:04:43.46ID:05wxJPem
RやSPSSって名前の試験まであるのね

490132人目の素数さん2018/08/22(水) 14:02:15.21ID:rznk0lAS
規約分数にするパッケージが探せなかったので自分でスクリプトを組んでみた。
エラー処置は省略。


reduce_fraction <- function(x,y){
a=x
b=y
r = a%%b # a=bq+r -> a%%b=b%%r
while(r){
a = b
b = r
r = a%%b
}
gcd=b
cat(x/gcd,'/',y/gcd,'\n')
invisible(c(x/gcd,y/gcd))
}


> reduce_fraction(2860,1082900)
11 / 4165

491132人目の素数さん2018/08/23(木) 01:12:39.77ID:0Wz4IoKN
rvest?でログイン必要なとこをスクレイピングする時って
login_page <- html_session("https://xxxxx")

login_form <- html_form(login_page)[[1]] %>%
set_values(AAA="xxxx", BBB="xxxx")

というので行けるはずなんですが
最終行のAAAとBBBってそれぞれhtmlタグのname=""のとこからとってくるんですよね?
これにここがname="user[email]"みたいに[]という記号はいってたらどうすればいいでしょうか?

492132人目の素数さん2018/08/25(土) 17:48:59.21ID:MdsDkupV
>>491
よろしくおねがいいしまーーーーーーーす

493132人目の素数さん2018/08/26(日) 22:01:57.06ID:2bj/HrEX
>>492
質問の意味がいまいち分からないから、誰も助言できないのでは?

> x <- 'user[email]'
> x
[1] "user[email]"

rvestパッケージを使ったことがないけど、
記号が入っていたらなぜ問題があるのかよく分からない。

494132人目の素数さん2018/08/26(日) 22:53:43.10ID:WnqFRbi9
というかどこのページのスクレイピングしたいか晒して
仮でもいいから

495132人目の素数さん2018/08/27(月) 06:16:23.20ID:x8E4FG0O
>>493
ちがいます
AAAがどれかよくみてみてください
””はありません

496132人目の素数さん2018/08/27(月) 06:19:43.00ID:x8E4FG0O
ちょっと直接はることはできないのでタグだけはります。
<input class="form-control" autofocus="autofocus" placeholder="Email address" required="required" type="text" value="" name="user[email]" id="user_email">>>494

497132人目の素数さん2018/08/27(月) 13:25:39.67ID:SdOxff6m
>>495
もっと分かるように説明しないと、追試できる情報も提供しないし。
もしかして変数名などに[]が入っている場合 にどうしたらよいかってこと?

> `user[email]` <- 1:5
> `user[email]`
[1] 1 2 3 4 5

498132人目の素数さん2018/08/27(月) 15:43:37.46ID:x8E4FG0O
>>497
多分行けました
超感謝!!

499132人目の素数さん2018/08/27(月) 20:20:32.13ID:5a+trmgU
ある問題のシミュレーションしようと思って
問題文の記号のまま

q <- function(x)
....

とやって
q(100)と入力すると
Rが終了することに気づいた。

500132人目の素数さん2018/08/27(月) 22:56:54.68ID:a+8R0Tu1
base::q()
を実行すると爆弾出るね

501132人目の素数さん2018/08/29(水) 22:45:10.53ID:BcwFyR33
ちょっとした疑問です。

空ベクトルの検出って長さ0以外に検出方法ってあるでしょうか?

> x=c(1,2)
> x=x[-1]
> x=x[-1]
> length (x)
[1] 0
> x
numeric(0)
> x==numeric (0)
logical(0)
> is.null(x)
[1] FALSE
> is.na(x)
logical(0)
> x==NULL
logical(0)
> x==NA
logical(0)
> length(x)==0
[1] TRUE

502132人目の素数さん2018/08/31(金) 22:46:15.16ID:IWQvY6FL
4点の座標を入力するとそれらを結ぶ四面体の体積を求めるスクリプトを書いてみた。
高さはパッケージnleqslvを使った近似計算。

# Calculate tetrahedron volume from cordinates
library(nleqslv)
Tetra <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
fn <- function(x,O,A,B,C){
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO # H on triangle ABC
AB=B-A
AC=C-A
c(HO%*%AB,HO%*%AC) # HO vertial to AB and AC
}
fn1 <- function(x) fn(x,O,A,B,C)
x=nleqslv::nleqslv(c(1/3,1/3),fn1)$'x'
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO
h=sqrt(sum(HO^2))
a=sqrt(sum((B-C)^2))
b=sqrt(sum((C-A)^2))
c=sqrt(sum((A-B)^2))
s=(a+b+c)/2
base=sqrt(s*(s-a)*(s-b)*(s-c))
V=1/3*base*h
return(V)
}

初期値は辺の長さ1の正四面体

> options(digits = 16)
> Tetra()
[1] 0.1178511301977579
> sqrt(2)/12
[1] 0.1178511301977579

多分、正常に動作していると思う。

503132人目の素数さん2018/09/01(土) 00:25:33.05ID:52Ub52jp
>>502
行列式det使うと簡単
> po <- c(1/2, sqrt(3)/6, sqrt(2/3))
> pa <- c(0,0,0)
> pb <- c(1,0,0)
> pc <- c(cos(pi/3), sin(pi/3), 0)
> det(cbind(pa-po,pb-po,pc-po))/6
[1] -0.1178511

504132人目の素数さん2018/09/01(土) 01:41:47.56ID:qG52f2Ee
ベクトルの三重積を教わったので、パッケージ pracma の外積crossを使った

tetrahedron <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
AO=A-O
BO=B-O
CO=C-O
as.numeric(abs(pracma::cross(AO,BO) %*% CO)/6)
}

4行で済んだ。

>>503

ありがとうございました。

パッケージに頼らずに計算できたのですね。

505132人目の素数さん2018/09/03(月) 18:55:17.77ID:S47YTHgP
データを解析する前にさらっと特徴を見たい時、皆さんはどんなコマンドを使っていますか?

私が思いつくのは
summary
boxplot
hist
pairs

です。こんなのも良いよってのがあったら教えてくださいm(_ _)m

※ライブラリの使用有無は問いません

506132人目の素数さん2018/09/03(月) 19:38:36.45ID:fowZfPON
>>505
BESTパッケージのplotPost
histに加えて95%CIを表示してくれる。

507132人目の素数さん2018/09/03(月) 20:05:14.66ID:OGj8hrn2
>>505
出力をデータフレーム型にしたいときはskimrパッケージ

508132人目の素数さん2018/09/03(月) 20:07:16.33ID:OGj8hrn2
あと、まだ試してないけsummarytoolsパッケージも面白そう

509132人目の素数さん2018/09/09(日) 22:35:22.23ID:9XY+z1xx
>>501
良い方法が見つからなかったので !length(x)で空白ベクトル判定とした。

文字列を逆に並べる再帰呼び出しスクリプト

> reverse <- function(x){
+ if(!length(x)) return(NULL)
+ c(Recall(x[-1]),x[1])
+ }
> cat(reverse(LETTERS[1:26]))
Z Y X W V U T S R Q P O N M L K J I H G F E D C B A

510132人目の素数さん2018/09/09(日) 23:21:54.06ID:1udV8jVZ
>>509
なぜ if (length(x) == 0) と書かないのか?
可読性って知ってるか?

511132人目の素数さん2018/09/10(月) 16:54:23.64ID:89eIPezd
>>510
文脈読めば
!length(x)で空白ベクトル
の流れ

512132人目の素数さん2018/09/10(月) 18:22:05.74ID:gLTPxqtw
>>511
length は数値を返す関数なのに、なぜわざわざ論理演算をするんだ?
プログラムは、ここの議論を知らない人が読む可能性のほうが高いんだぞ。

513132人目の素数さん2018/09/10(月) 19:12:52.95ID:ZYY4OYkH
>>512
正規分布が一様分布より大きい期待値の算出に
mean(rnorm(100)>runif(100))
って書かない?
俺は
sum(ifelse(rnorm(100)>runif(100),1,0))/100
って書いたりしないけど。

514132人目の素数さん2018/09/10(月) 19:43:25.60ID:SKQ8XoAj
>>512
え、常套手段やん?!

515132人目の素数さん2018/09/10(月) 19:50:35.36ID:ZYY4OYkH
>>514
論理値を数値に置き換えて計算しているから
数値を論理値にしても別に違和感がない。
可読性は慣れの問題。

516132人目の素数さん2018/09/10(月) 19:54:17.32ID:SKQ8XoAj
>>515
アンカー間違えとる?
だから、そんなの常套手段やん?

517132人目の素数さん2018/09/10(月) 21:30:23.37ID:gLTPxqtw
>>513
それとこれとは話が別だ。
logical は TRUE か FALSE の二値しかとらないから、それから数値への変換は自明。
多種の値をとる数値に論理演算をするのはいただけない。

518132人目の素数さん2018/09/10(月) 21:36:20.92ID:gLTPxqtw
>>515
is.empty.vector のような関数があるならそれでよいが、length を使っているのだから、それを勝手にベクトルが空かどうかの判断に使っているのはあなたであって、他人はそのようには考えない、といっているのだ。
自分だけしかこーどを見ないならべつに構わないが、こういうところに晒すのはよくない。
ましてや常套手段などと言って他人に教え込むのはやめてもらいたい。

519132人目の素数さん2018/09/10(月) 21:51:19.40ID:5QS5/GHY
>>518
常套手段というのはmeanの話じゃね?

520132人目の素数さん2018/09/10(月) 21:55:50.32ID:5QS5/GHY
whileの中は1でも2でもよくね?

n=0
while(1){
if(n>10) return(10)
n=n+1
}


n=0
while(2){
if(n>10) return(10)
n=n+1
}

521132人目の素数さん2018/09/10(月) 22:00:57.10ID:5QS5/GHY
>509は
配列を逆順に並べる再帰呼出しのコード。

Rにはrevという関数がある。

そのコードみてみ!

lengthを真偽判定に使ってますがな。


> base:::rev.default
function (x)
if (length(x)) x[length(x):1L] else x

522132人目の素数さん2018/09/10(月) 22:18:05.27ID:DJR6mPp+
>>518
>521みた?
if(length(x)!=0)になってないよ。

523132人目の素数さん2018/09/11(火) 00:19:39.45ID:WPUC3B84
>>518
>こういうところ
誰も鵜呑みにしない便所の落書き。。。

524132人目の素数さん2018/09/11(火) 08:48:42.34ID:QUqp/jpE
!を引数が数値のときは0か否かを返す関数と理解すれば
if(!0) print(!1)の結果もサクッとわかる。

525132人目の素数さん2018/09/11(火) 09:08:32.20ID:QUqp/jpE
数値を非零かどうか論理値に変換して処理するかは
関数によるな。

if やwhileは変換処理しているけど
whichだとエラーになるな。

> which(!0)
[1] 1
> which(0)
Error in which(0) : argument to 'which' is not logical
> which(2)
Error in which(2) : argument to 'which' is not logical
> which(!3)
integer(0)
>

526132人目の素数さん2018/09/11(火) 09:12:56.04ID:QUqp/jpE
!ってベクトル対応しているな
> which(!(-1:1))
[1] 2
> !(-2:2)
[1] FALSE FALSE TRUE FALSE FALSE

527132人目の素数さん2018/09/15(土) 08:05:49.52ID:Vl7XZ52q
!をnotではなくis.falseとかis.zero変換すれば( ・∀・)イイ!!

528132人目の素数さん2018/09/15(土) 16:42:00.50ID:h0gUCZ3r
>>527
is.zero だと?
だったら最初から == 0 てよいじゃないか。

529132人目の素数さん2018/09/15(土) 16:45:38.59ID:h0gUCZ3r
>>525
だいたい、if や while は関数じゃないし。
while (1) などと書くのはCに毒されているんじゃねーの?
Rなら while (TRUE) のほうが良いし、repeat というのもある。

530132人目の素数さん2018/09/15(土) 16:45:41.31ID:Vl7XZ52q
>>528
スクリプトを読むときに脳内変換する癖をつけておけば
可読性が悪いと思わずにすむ。

531132人目の素数さん2018/09/15(土) 16:46:13.42ID:Vl7XZ52q
>>518
>521みた?
if(length(x)!=0)になってないよ。

532132人目の素数さん2018/09/15(土) 17:19:35.45ID:h0gUCZ3r
>>531
Rのライブラリ?の書き方なんてあてにならないよ。
そもそも、rev も 509 の reverse も演算の回数を必要最低限にするなら、条件は length(x)) >= 2 または > 1 だ。
まあ if (length(x)) のほうが速かったのかも知れないが、常にそうなるとは限らないのだから、早すぎる最適化は諸悪の根源だ。

533132人目の素数さん2018/09/15(土) 19:30:34.55ID:Vl7XZ52q
可読性の次は速度かよ。
自分の流儀と違う { の使い方でも文句いいそう。

534132人目の素数さん2018/09/15(土) 19:44:27.43ID:VibLIqgl
0,1を1000万個作って
!
!=
>
==
での真偽判定を各々10回する時間を出してみた。


> x=rbinom(1e7,1,0.5)
> system.time(replicate(10,!x))
user system elapsed
1.27 0.33 1.59
> system.time(replicate(10,x!=0))
user system elapsed
2.39 0.36 2.75
> system.time(replicate(10,x>0))
user system elapsed
2.58 0.31 2.92
> system.time(replicate(10,x==1))
user system elapsed
2.60 0.36 2.99

535132人目の素数さん2018/09/15(土) 19:57:45.72ID:VibLIqgl
!x と x!=0 で10回やったが、

> f <- function(){
+ x=rbinom(1e7,1,0.5)
+ (a=system.time(!x)[3])
+ (b=system.time(x!=0)[3])
+ a<b
+ }
> (re=replicate(10,f()))
elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
elapsed
TRUE

結果は再現された。

536132人目の素数さん2018/09/15(土) 20:06:37.73ID:VibLIqgl
lengthが絡むと

> x=1:1e8
> f1 = function(x) if(!length(x)) return(NULL) else x=x[-1]
> f2 = function(x) if(length(x)==0) return(NULL) else x=x[-1]
> system.time(f1(x))
user system elapsed
2.38 0.53 2.93
> system.time(f2(x))
user system elapsed
2.05 0.61 2.69

御指摘の通り、!の負け

537132人目の素数さん2018/10/04(木) 17:02:48.04ID:247Ted9r
>>444
>遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
>詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ

Ubuntu 18.04.1 LTSにアップグレードしたけどできない
いままでスクレイピングにつかってrvestのログインコマンドもエラー出るようになったし最悪だ。。。。

538132人目の素数さん2018/10/04(木) 20:04:49.68ID:W6fWyA5T
アップグレードってことは文字入力がFcitxのままになってない?Fcitxならパッチ要だよ。
Voyager 18.04LTS(xubuntu 18.04LTS)ではパッチなしでRStudioに日本語入力できた。

539132人目の素数さん2018/10/26(金) 12:58:59.43ID:8LsI2NoU
昨日、PCを新しくして最新のR入れたんだけど、
作業ディレクトリの固定化できなくね?
前まではショートカットのプロパティで作業フォルダ指定してたんだけど、
新しいR、何度やっても作業ディレクトリがデフォルトのドキュメントに戻りやがる

どうすればいいか、誰がご教示を。

540132人目の素数さん2018/10/26(金) 13:04:50.23ID:I9GOHEF6
>>539
プロパティでコマンドラインオプションを消せばよい。
--cd-ほにゃらら
というのが付いてるでしょ。

541132人目の素数さん2018/10/26(金) 13:13:30.11ID:8LsI2NoU
>>540
あ、なるほど
ありがとうございます、多謝

542132人目の素数さん2018/10/26(金) 23:05:34.11ID:5yNEfvRx
xtsの取り扱い方について質問失礼します。
以下のコードで、dataには2018/10/20 00:00〜2018/10/25 00:00の1時間刻みのaとbの値のxts型オブジェクトが入りますが、(aとbは時系列データのイメージです)
このdataの毎日09:00の値を添え字を使って取り出す方法はありますか?
(data["2018-10-21 00:00"]とすると10/21の00:00が取り出せるような感じで)

library(xts)
time <- seq(as.POSIXct("2018-10-20 00:00"), as.POSIXct("2018-10-25 00:00"), by="1 hour")
a <- rnorm(length(time))
b <- rnorm(length(time))
data <- xts(x = cbind(a,b), order.by = time)

思いついたのはfor文を使って以下のようにするのかなと思ったのですがもう少しいい方法があるだろうなと思いまして…

data_result <- data[10] # まずxts型の入れ物を作る(2018/10/20 09:00)
for (i in 1:4) {
data_result <- rbind(temp,data[10+24*i]) # 2日目以降はこれにくっつける
}
data_result

543132人目の素数さん2018/10/26(金) 23:12:48.52ID:Xg1/ChR5
失礼、>>542の下のコードはこっちです

data_result <- data[10] # まずxts型の入れ物を作る(2018/10/20 09:00)
for (i in 1:4) {
data_result <- rbind(data_result, ,data[10+24*i]) # 2日目以降はこれにくっつける
}
data_result

544132人目の素数さん2018/10/26(金) 23:38:29.60ID:oEWqNGH8
>>543
idx <- seq(10, length(time), by=24)
result <- data[idx]
じゃダメなのか?

545132人目の素数さん2018/10/27(土) 01:10:34.48ID:NQVW6zPX
>>544
ああ、これがいいですね!お恥ずかしい…
ありがとうございます

546132人目の素数さん2018/11/01(木) 22:44:27.08ID:xCdOvDq8
巨大数を扱えるというふれこみのRmpfrって正確じゃないな。
50の階乗を計算させてみた

R with Rmpfr
> mpfr(factorial(50),1e5)
1 'mpfr' number of precision 100000 bits
30414093201713018969967457666435945132957882063457991132016803840

Haskell:
Prelude> product[1..50]
30414093201713378043612608166064768844377641568960512000000000000

Python:
import math
print(math.factorial(50))
30414093201713378043612608166064768844377641568960512000000000000

Wolfram:
https://www.wolframalpha.com/input/?i=50!
30414093201713378043612608166064768844377641568960512000000000000

547132人目の素数さん2018/11/02(金) 23:29:01.53ID:fMme23mF
こまけぇこたぁいいんだよ!!

548132人目の素数さん2018/11/03(土) 09:32:49.16ID:NiwdVATo
R3.5.xから、Windows環境だと日本語パスまわりでエラーがでてどうにもならない
enc2native()しても駄目だし
とりあえずcairo_png()だけでも…

549132人目の素数さん2018/11/03(土) 10:31:12.42ID:pqZePX+I
一度つかって日本語周りの処理でエラー出てきたんで未だに3.4.4浸かってる

550132人目の素数さん2018/11/04(日) 21:55:16.95ID:lTCeMsqQ
統計の質問はこのスレでいいんだろうか?

ある医院に1時間あたり平均5人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
1時間あたりの平均診療人数は6人で、一人あたりの診療時間は指数分布に従うとする。
診察までの平均の待ち時間は何時間か?

このシミュレーション解はこれであってる?

N=1e6
λ=5
μ=6
sum(rpois(N,λ)*rexp(N,μ))/N

> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.8331036

60かけて、分にすると
[1] 49.86565

551132人目の素数さん2018/11/08(木) 12:35:54.73ID:wKTjJ6Fa
>>306
grep 使って書いてみた。
# mhs(c(1,0,1,1,0,1,1,1)) : return 3
mhs = function(x){ # maximum head sequence
y=paste(x,collapse='')
str="1"
if(!grepl(str,y)) return(0)
else{
while(grepl(str,y)){
str=paste0(str,"1")
}
return(nchar(str)-1)
}
}

system.time(mean(replicate(1e4,mhs(sample(0:1,100,rep=T))>=5)))


# N(=100)回コインをn(=5回)以上続けて表がでるか?TRUE or FALSE
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}

system.time(mean(replicate(1e4,seqn())))


結果は、
> system.time(mean(replicate(1e4,mhs(sample(0:1,100,rep=T))>=5)))
user system elapsed
4.56 0.01 4.61
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
1.05 0.00 1.07

逆に4倍時間がかかるようになった

552132人目の素数さん2018/11/08(木) 14:05:00.34ID:1h1GfW+d
今さらだけどそれってrle使っちゃだめなの?

553132人目の素数さん2018/11/08(木) 15:13:38.19ID:UNhl34kg
coinの表裏を1と0で表して、その累積和を取ったベクトルを用意して
そのベクトルの5個前とのdiffの要素の中に5が一回でも現れることが、一回でも5回連続で表が出たことがあることの必要十分条件
grepはどうしても時間がかかるしif文もできれば使いたくない

searchseq <- function(n=5,N=100,p=0.5,trial=10000){
result <- 0 # 表が5回以上出た回数を数える入れ物
for (i in 1:trial){
coin <- rbinom(N,1,p)
coincumsum <- cumsum(coin)
coindiff <- diff(coincumsum,5)
#diff(coincumsum,5)の要素に一個でも5があれば表が5回以上出たことがあったということ
#anyでT/Fにして、sumで0/1にする
result <- result+sum(any(coindiff==5))
}
return(result/trial)
}

結果もよさそう
> searchseq()
[1] 0.80793

> system.time(mean(replicate(10000,mhs(sample(0:1,100,rep=T))>=5)))
ユーザ システム 経過
0.92 0.00 0.93
> system.time(mean(replicate(10000,seqn())))
ユーザ システム 経過
0.29 0.00 0.30
> system.time(searchseq())
ユーザ システム 経過
0.21 0.00 0.20

554132人目の素数さん2018/11/08(木) 22:10:55.78ID:I6htNxzg
あっ、コード中の5はnに置き換えてね

555132人目の素数さん2018/11/08(木) 23:32:48.68ID:/ZlhhGjJ
>>553
レスありがとうございました。

rleとも比べてみました。

> system.time(mean(replicate(1e4,seqn())))
user system elapsed
4.290 0.000 4.377
>system.time(mean(replicate(1e4,max(rle(rbinom(100,1,0.5))$len)>=5)))
user system elapsed
3.620 0.000 3.742
> system.time(mean(replicate(1e4,mhs(rbinom(100,1,0.5)>=5))))
user system elapsed
2.390 0.000 2.435
> system.time(searchseq())
user system elapsed
1.880 0.000 1.988

556132人目の素数さん2018/11/08(木) 23:35:10.38ID:/ZlhhGjJ
>>555
greplの逆転はsampleをrbinomに変えたためでしょう。

> system.time(replicate(1e5,sample(0:1,100,rep=T)))
user system elapsed
7.200 0.180 7.591
> system.time(replicate(1e5,rbinom(100,1,0.5)))
user system elapsed
5.980 0.230 6.319

557132人目の素数さん2018/11/09(金) 00:20:12.33ID:BZRFS9lT
不勉強ながらrle関数って初めて知ったけど使いやすそうだな

558132人目の素数さん2018/11/09(金) 00:36:46.75ID:Ui6BpICy
>>555
rleは0も数えているから間違っているんじゃ?

559132人目の素数さん2018/11/09(金) 02:01:36.95ID:Qla0VxTD
>>558
ご指摘ありがとうございました。
修正しました。
rle1 <- function(N=100,n=5,p=0.5){
r=rle(rbinom(N,1,p))
max(r$len[which(r$val==1)])
}
結果
> system.time(mean(replicate(1e4,max(rle1()>=5))))
user system elapsed
4.430 0.000 4.546
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
4.490 0.010 4.609
> system.time(mean(replicate(1e4,mhs(rbinom(100,1,0.5))>=5)))
user system elapsed
7.440 0.000 7.656
> system.time(searchseq())
user system elapsed
1.950 0.000 2.066
>

560132人目の素数さん2018/11/09(金) 07:25:00.97ID:Qla0VxTD
無理矢理1行にして実行

system.time(mean(replicate(1e4,any(diff(cumsum(rbinom(100,1,0.5)),5)==5))))
user system elapsed
1.820 0.000 1.886
>
> system.time(mean(replicate(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[wh
<e(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[whi ch(values==1)])>=5))))
user system elapsed
4.370 0.010 4.478

561132人目の素数さん2018/11/09(金) 07:59:47.47ID:rijPDFSR
>>560
意味もなくforループ回してた上に毎回sum使って真偽値を数値に変換してたけど
replicate使って最後に一回だけmean取ると2.066→1.886で1割短くなるのね
他人のコード読むのは勉強になる

562132人目の素数さん2018/11/09(金) 09:06:56.99ID:Ui6BpICy
>>559
whichいらねーよ。

563132人目の素数さん2018/11/09(金) 09:44:07.41ID:5AnUTlVm
>>559
全部が0のとき、エラーになるので修正

rle01 <- function(x){ # c(0,1,1,1,0,0) => return 3
if(sum(x)==0) return(0) #c(0,0,0,0,0,0) => return 0
else{
r=rle(x) # Run Length Encoding
max(r$lengths[which(r$values==1)]) # max length of value 1
}
}

動作確認

> rle01(x<-rbinom(100,1,0.5)) ; x
[1] 8
[1] 1 1 0 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1 0 0 0 1 1 1 1 0 0 1 1 1 1 1 0 0
[36] 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 0 1 1 0 1 0 1 0 0
[71] 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 0 0 1 0 0 1 1 1 0 1 1
> rle01(x<-rbinom(100,1,0)) ; x
[1] 0
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

564132人目の素数さん2018/11/09(金) 10:10:01.39ID:Ui6BpICy
>>563
センスないなあ。
rle01 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}

565132人目の素数さん2018/11/09(金) 10:40:31.84ID:5AnUTlVm
>>564
whichは余分でした。

sumの方がrleより高速だと思ったらから、すべて0の場合はrleを呼ばないことにしただけ。
0連続でやるとこれだけ差がつく。

rle01 <- function(x){ # c(0,1,1,1,0,0) => return 3
if(sum(x)==0) return(0) #c(0,0,0,0,0,0) => return 0
else{
r=rle(x) # Run Length Encoding
max(r$lengths[which(r$values==1])] # max length of value 1
}
}
rle012 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}

> x=rep(0,1e8)
> system.time(rle01(x))
user system elapsed
0.3 0.0 0.3
> system.time(rle012(x))
user system elapsed
7.36 4.52 13.25
>

566132人目の素数さん2018/11/09(金) 10:52:23.39ID:Ui6BpICy
>>565
そんな特別な場合のことに対して高速化するのは愚の骨頂。
1が一つでもある場合にsumを呼ぶのは余計じゃないのか?
余計なwhichがあるくらいだから、まずは素直にやりたいこと、やるべきことを正しく書くようにしたら?

567132人目の素数さん2018/11/09(金) 11:01:40.74ID:5AnUTlVm
1000本に1本あたる宝クジを100本買って続けて2本あたる確率のシミュレーション解の算出時間比較。
確率の理論値は9.8897353347449091e-05

> system.time(mean(replicate(1e4,rle01(rbinom(N,1,p))>=n)))
user system elapsed
0.61 0.00 0.64

> system.time(mean(replicate(1e4,rle12(rbinom(N,1,p))>=n)))
user system elapsed
1.97 0.00 2.03

568132人目の素数さん2018/11/09(金) 15:22:38.24ID:5AnUTlVm
>>567
分数で表すと

890788167367/9007199254740992

9.889735334744909e-05

569132人目の素数さん2018/11/09(金) 19:43:44.34ID:0M0agNOf
JPXのデータ、ファイル形式csvを読み込もうとするとうまく行かないんですが
どんな引数をつければいいですか

570132人目の素数さん2018/11/10(土) 10:36:49.41ID:QJ6NJqU7
>>566
実はシミュレーションじゃなくて
漸化式からのプログラム解を分数表示するプログラムはpythonで作成済。
ここに置いた。
https://egg.2ch.net/test/read.cgi/hosp/1540905566/77

571132人目の素数さん2018/11/10(土) 11:20:41.80ID:QJ6NJqU7
コインを100回ふったときの表連続の最大数が5であったときの
このコインの表がでる確率の期待値、モード比、信頼区間を求めるのが次のネタ。

unirootで算出できたけど
シミュレーションはどうすればいいのかアイデアが浮かばない。
MCMCで解決できるかなぁ?

これをシミュレーションで検証したい。

$Rscript main.r
lower mean mode upper
0.2487456 0.4469764 0.4589692 0.6386493

http://tpcg.io/asKRE9

572132人目の素数さん2018/11/11(日) 20:11:21.81ID:ODPKEEGK
1億回のコイントスで何回連続して表がでる確率が高いかRでやってみた。

# maximal sequential head probability at 10^8 coin flip
> y
[1] 2.2204460492503131e-16 2.2204460492503131e-16
[3] 8.8817841970012523e-16 1.5543122344752192e-15
[5] 3.5527136788005009e-15 6.8833827526759706e-15
[7] 1.4210854715202004e-14 2.8199664825478976e-14
[9] 5.6843418860808015e-14 1.1346479311669100e-13
[11] 2.2737367544323206e-13 4.5452530628153909e-13
[13] 9.0949470177292824e-13 1.8187673589409314e-12
[15] 3.6379788070917130e-12 7.2757355695785009e-12
[17] 1.4551915228366852e-11 2.9103608412128779e-11
[19] 5.8207660913467407e-11 1.1641509978232989e-10
[21] 6.6493359939245877e-06 2.5720460687386204e-03
[23] 4.8202266324911647e-02 1.7456547460031679e-01
[25] 2.4936031630254019e-01 2.1428293501123058e-01
[27] 1.4106434838399229e-01 8.1018980443629832e-02
[29] 4.3428433624081136e-02 2.2484450838189007e-02

25回が続くのが4回に1回あることになる。
pythonで25回以上と25回ちょうどになるのを計算させてみた。

その結果、
Over 25
6977459029519597/9007199254740992
= 0.7746535667951356
Just 25
2246038053550679/9007199254740992
= 0.24936031612362342

高速化を狙ってCに移植したら
100万回で暴走。
https://egg.5ch.net/test/read.cgi/hosp/1540905566/132

573132人目の素数さん2018/11/12(月) 21:04:57.19ID:boi8bvdM
"マラソン大会の選手に1から順番に番号の書かれたゼッケンをつける。
用意されたゼッケンN(=100)枚以下の参加であった。
無作為に抽出したM(=5)人のゼッケン番号の最大値はMmax(=60)であった。
参加人数推定値の期待値とその95%信頼区間を求めよ"

decken <- function(M, Mmax, N, conf.level=0.95){
if(Mmax < M) return(0)
n=Mmax:N
pmf=choose(Mmax-1,M-1)/choose(n,M)
pdf=pmf/sum (pmf)
mean=sum(n*pdf)
upr=n[which(cumsum(pdf)>conf.level)[1]]
lwr=Mmax
c(lower=lwr,mean=mean,upper=upr)
}
decken(M=5,Mmax=60,N=100)

> decken(M=5,Mmax=60,N=100)
lower mean upper
60.0000 71.4885 93.0000

これをシミュレーションで確認したい。

# simulation
M=5 ; Mmax=60 ; N=100
sub <- function(M,Mmax,N){
n=sample(Mmax:N,1) # n : 参加人数n
m.max=max(sample(n,M)) # m.max : n人からM人選んだ最大番号
if(m.max==Mmax) return(n)
}
sim <- function(){
n=sub(M,Mmax,N)
while(is.null(n)){
n=sub(M,Mmax,N) # 最大番号が一致するまで繰り返す
}
return(n)
}
runner=replicate(1e4,sim())
summary(runner) ; hist(runner,freq=F,col="lightblue")
quantile(runner,prob=c(0,0.95))
cat(names(table(runner)[which.max(table(runner))]))

> summary(runner) ; hist(runner,freq=F,col="lightblue")
Min. 1st Qu. Median Mean 3rd Qu. Max.
60.00 63.00 68.00 71.43 77.00 100.00
> quantile(runner,prob=c(0,0.95))
0% 95%
60 93
> cat(names(table(runner)[which.max(table(runner))])) # 最頻値
60

結果は確認できたけど、もっと高速なシミュレーションアルゴリズムはあるだろうか?

新着レスの表示
レスを投稿する