【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)
}

新着レスの表示
レスを投稿する