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

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言語】統計解析フリーソフトＲ 第５章【GNU R】
http://rio2016.2ch.net/test/read.cgi/math/1380168442/

板違う

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

耳栓をしたら世界が変わってワロタ

￥

耳栓をしたら世界が変わってワロタ

￥

耳栓をしたら世界が変わってワロタ

￥

耳栓をしたら世界が変わってワロタ

￥

耳栓をしたら世界が変わってワロタ

アンドロイドにRをインストールしてみた。
テキストベースなら簡易統計電卓として使える。
依拠したのはこれ。
https://stackoverflow.com/questions/36968411/installing-r-on-android
これくらいの簡単な計算ならタブレットでできるようになった。
いちいちパソコンを立ち上げなくてもすむ。

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

>>297
クリップボード経由でスクリプトを実行してみたが
日本語が混ざるとエラーになるな。
フィーリングカップル5対5のカップル数とカップル成立確率の
シミュレーションを走らせてみた。

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

きちんと動く。

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.
ということでパッケージのインストールはエラーになった。

>>25
>dependency hell
具体性がない。どこ情報？

>>26
http://www.r-ohjelmointi.org/?p=1434

Console R はコンパイラーが使えたからPremium版を買ったけど
今は公開中止で入手できない。

>>27
Debianの話ではなくて、Andoridか、紛らわしい。
ってか、クロスコンパイル環境でビルドするのではなくて、
DebianのバイナリをAndoridの中のエミュレータで動かすという話？

こいつらスマホからR使ってんの？

FireHDで使いたい

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

>>31
先ず、Google Playをいれて>23の手順で使えるようになるはず。

>>29
GNURoot Debian provides a method for you to install and use Debian and its associated applications/packages alongside Android.
とあるからエミュレータなのだろう。

>>31
他の泥タブレットと同じくFireHDでもできたよ。

気づいたら前スレ埋められてたのか
いつも思うがRに親でも殺されたのかと

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

2^127-1 =170141183460469231731687303715884105727は素数

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

Ｒだと

> 2^127-1
[1] 170141183460469231738248242066868668888

と表示された。

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

2^1279-1も素数

10407932194664399081925240327364085538615262247266704805319112350403608059673360
29801223944173232418484242161395428100779138356624832346490813990660567732076292
41295093892203457731833496615835504729594205476898112116936771475484788669625013
84438260291732348885311160828538416585028255604666224831890918801847068222203140
521026698435488732958028878050869736186900714720710555703168729087

なんかパッケージないの？

gmpというパッケージがあったので使ってみた。

メルセンヌ素数算出のスクリプトを組んでみた
2^n-1 で　ｎ　<　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

Rって何のための言語だっけ

Sのフリー版

Rに親を殺されたマンはよ

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

>>60
貧乏人のためのデータプログラミング言語

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

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


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

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


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


☆数学コテ紹介☆

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

耳栓をしたら世界が変わってワロタ

￥

￥

￥

￥

￥

￥

￥

￥

耳栓をしたら世界が変わってワロタ

￥

耳栓をしたら世界が変わってワロタ

￥

耳栓をしたら世界が変わってワロタ

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

耳栓をしたら世界が変わってワロタ

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

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

>>356
実行すると
> LE(f,96)
[1] 3.4
96歳女性の平均余命は3年以上あるらしい。
死因は多臓器不全にする方が正しかったのだろうか？

> LE(m,100)
[1] 2.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を参照。

どうでもいいけど、平成27年は完全生命表の年なのに、
なぜわざわざ簡易生命表を使うの？

>>111
5〜6年ほど前にやってたな。

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

>>113
公表されていればurl希望

>>115
第22回生命表（完全生命表）
ttp://www.mhlw.go.jp/toukei/saikin/hw/life/22th/index.html

なお、最新版の平成28年の簡易生命表はこちら
ttp://www.mhlw.go.jp/toukei/saikin/hw/life/life16/index.html

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

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

>>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

微妙に結果が違うな。

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

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)

> lottery(100,1)
2.5% 5% 50% 95% 97.5%
16 22 70 97 98
mean
66.31167

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)

Linux版のRでもWindows版みたいに「入力:青字、出力:赤字」みたいなインターフェースにする設定ってある？

>>125
coloroutパッケージ入れて、色を割り当てるとできるかも。
アーカイブに入っちゃっているけど、使える。
setOutputColors256()のヘルプを参照。

version 3.4.2になってた。

順番を表示させる、*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"

>>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"

11,12,13だけ除外かな

>>128
なんでそこで && を使うかなあ。&& の意味わかってる？
ifelse には & を使うことになるのが普通だし、そうでなければ if 〜 else のほうがよい。

>>129
ご指摘とパッケージの御教示ありがとうございました。

>>131
?&で勉強になりました。
御教示ありがとうございました。

>>130
toOrdinal のコードをみたら　11,12,13は別扱いしてあった。

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

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万で合致を確認。

いろいろとご助言ありがとうございました。　<(_ _)>ﾍﾟｺﾘ

パチンコのスレにこういう記載があった。

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

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

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

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

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

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

分散が変わる？

シミュレーションすると

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

お知恵を拝借したい。

よくわからん。rbinomとか使えば？

>>139
俺もそう思った

自分なら、こんな感じ
> 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)
同じ結果だよ

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)
結果は同じ

>139-142

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

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

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

>>142
大したことではないけど

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

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

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

合計の期待値の標準偏差はAの本数が増えると減るね

宝くじを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

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

>>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

別の話題提供。

　ゴルゴ13は100発100中

　ゴルゴ14は10発10中

　ゴルゴ15は1発1中

とする。

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

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

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

length(b2)

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

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

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

どうでもいいことだけど、

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

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

>>126
125です。レス遅くなったけど、やりたいことができました。ありがとう。

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なんです。

>>154
>3行2列
さっぱりわからん。
例えば、結果の[1,1,1]の5は、どれとどれの平均？

>>155
6行4列の左上の6つ1,2,3,7,8,9の平均5だと思う。

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

>>155
>>154です。
説明が悪くてすみません。
>>156の通りです。

>>157
ありがとうございます。
手元にデータがないので明日試してみます。

なるほど、じゃあその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

あぁ、すでに >>157 が回答を書いていたorz

> 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
でよくね？

最初の1行で

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

エラーがでたｗ

正しいエラーｗはこっち

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

そもそも、10000x10000x400だと10000行は3行にならないな。最終行を無視して9999行にするのかな？

>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

>

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

先輩方，>>154です。
巨細にわたってご教示・ご議論いただきありがとうございました。
>>166のコードを参考にして，無事処理することができました。
実行速度も申し分なく，効率的に作業が進められそうです。
今回は私には目から鱗な解法で，勉強になりました。
質問して良かったと思っています。。本当にありがとうございました。

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

処理速度を比べてみた。

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氏のスキルに改めて感服。

# 既約剰余類群 (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

# 既約剰余類群での加減乗除
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

名前も一致して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

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)

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歳なので、老衰にせずに
主治医指定の病名で診断書を作成。
カルテ記載のしっかりした患者の診療は負担なくこなせて（・∀・）ｲｲ!!

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

>>174
そりゃあ、outer に渡す関数がまずいんだろ。
outer に渡す関数は、完全にベクトル対応になっていなければならないんじゃないの？
でなければ、わざわざ outer 使うこともなかろう。

>>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

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

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

>>176
outer から呼び出されたときの引数を見てみればわかるよ。たぶん1度しか呼ばれず、できあがりの行列の要素すべてを計算するための引数がベクトルで渡されているはずだから。

>>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, ...)

# 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

ユークリッドの互除法をつかうとこんな感じ

> 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

文系研究者のワイ、困惑
ディープな世界やでほんま

#.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)
}

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

こういう問題を考えるのは楽しいね。

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

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

ベクトルは再利用されるんだな


> 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

> 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

Brunner-Munzel検定って比較する両群に重なりがないとｐ値がＮＡになるんだが、
これってバグなのか？



> 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

突然ですが質問があります。最近Rソフトを利用し始めた者のです。単刀直入に、
ボロノイ図を制作した場合に分割された面積の平均を出力するには、どうすればよいのでしょうか？

現状としては住所をグーグルマップの緯度経度に変換し、分割まではできています。
宜しければどなたかお力添えをいただけないでしょうか？

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

>>198
SpatialPolygonsクラスまでできるといると仮定して、
各ポリゴンの中にareaスロットがあるから、
それを抽出して、その平均を取れば良い。
ボロノイ図がSpatialPolygonsクラスになっていないなら、
まずは変換。

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

ボロノイ図の各領域の面積の平均値って、単に
描画領域の総面積÷点数 じゃないの？
何か暗黙の前提があるの？

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

Stan って　離散量をparameter 扱いできないので

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

JAGSだと

m 〜　dcat(θ）

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

この理解でいいかな？

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

以前、ボロノイ図について質問させていただいた初心者です。
質問後色々と試行錯誤したのですが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
と入力しボロノイ図は出せているのですが、

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

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

>>242
検証用のデータがないと誰も追試しないと思う。
>154は模擬データがあったのでやる気になった。

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

関数のオプションの取り得る値を全て表示する方法って知らない？
例えばcor関数のmethodオプションの取り得る値は"pearson", "kendall", "spearman"みたいに。
そういうのは無いのかな？
初めて使う関数のオプションいちいちググるの面倒くさい。

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

>>265
ヘルプ表示じゃ駄目なのかい？

>>265
> args(cor)
function (x, y = NULL, use = "everything", method = c("pearson",
"kendall", "spearman"))
NULL

>>276
前にヘルプで見ようと思ったら書いてなかった覚えが…。マニアックな関数だったんで、そんなもんかな〜と諦めてた。やっぱ基本はヘルプか。
>>277
argsは見れる関数と見れない関数があるのであんまり役に立たない印象。

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

>>278
>argsは見れる関数と見れない関数があるのであんまり役に立たない
だったら、例が悪いよ。見られない関数の例をあげないと

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

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

rm(list=ls())

で　R　のメモリーをクリア

先週長々と質問させていただいたものです。その後、試行錯誤して行った結果十解決しました。

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

N(=100)回コインをなげてｎ(=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

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

ポアソンでやるやつだっけ

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()

>>303
単なる二項分布。
コインが５回続けて表がでたら、0.5^5 <0.05なのでイカサマコインといわれちゃいそうなんだが、
１００回やってみると案外、５回表が続くので確率を計算しようと思ったが、解析的にできる頭がないので
シミュレーションしてみた。

>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

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

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

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

すまん、教えてほしいだけど
分析するために初めてＲをインストールしようと思って、このスレのあるように公式サイト行ったら、esetが「JS/Redirector.NAV トロイの木馬」を検知したんだが…；
どうしたらいいだ…

そっちじゃなくてCRAN行け

ホームページダメとかわけわかんねぇ・・・とんでもねぇな
ありがとう

>>307
エロサイトにアクセスしてないｗ？

>>309
CRANはここ
https://cran.r-project.org


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

>>310
してねぇよｗｗ
totalvirusで調べると一件引っかかるし、Redirector検知だから、しょうがないね

col=rgb(runif(1),runif(1),runif(1),runif(1))で色指定すると

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

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

stanやJAGSのコードでgamma関数を使おうとして
y = gamma(x)
と、書いたらエラーになった。

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

ある大学の入学者男女の比率は１であるという帰無仮説を検定する課題が花子と太郎に課された。

花子は５０人を調査できたら終了として入学者を５０人をみつけて１８人が女子であるという結果を得た。
帰無仮説のもとで
５０人中１８人が女子である確率は 0.01603475
これ以下になるのは５０人中０〜１８人と３２〜５０人が女子の場合なので
両側検定して
> 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

一方、本番と十八番が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は0.005772512
これ以下になるのは２３人以下５０人以上番めで女子１８人めがみつかった場合なので
両側検定して
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
で帰無仮説は棄却される。

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

花子の横軸は裏口入学者数、太郎の横軸はサンプル数なので
サンプルでの裏口入学率を横軸にして９５％信頼区間を示す。
花子の検定での信頼区間は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
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。

（訂正）
ある大学の入学者男女の比率は１であるという帰無仮説を検定する課題が花子と太郎に課された。

花子は５０人を調査できたら終了として入学者を５０人をみつけて１８人が女子であるという結果を得た。
帰無仮説のもとで
５０人中１８人が女子である確率は 0.01603475
これ以下になるのは５０人中０〜１８人と３２〜５０人が女子の場合なので
両側検定して
> 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

一方、十八という数字が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は0.005772512
これ以下になるのは２３人以下５０人以上番めで女子１８人めがみつかった場合なので
両側検定して
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
で帰無仮説は棄却される。

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

花子の横軸は女子数、太郎の横軸はサンプル数なので
サンプルでの女子の割合を横軸にして９５％信頼区間を示す。
花子の検定での信頼区間は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
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。

kainokousiki<-function(a,b,c){return (-b+sqrt(b^2-4*a*c))/(2*a)} #解の公式
kainokousiki(1,-5,6)

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

自己解決
かっこが足りなかった

Pythonのスレはないのか

￥

￥

￥

￥

￥

￥

￥

￥

￥

￥

機械学習をきっかけにPythonに逆転された感じだね

『Rを使った〜』だとPythonじゃないのかよって思うよね

ずっとRだけでPython触ったこと無いけど、覚え直す価値ある？
環境構築からもう面倒なイメージ

Python自体は:と直後のインデントさえ気を付ければかなり簡単
3系は数が全て小数扱いなので楽

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

アナコンダてのがRぽくできるのね、ありがとう
dplyrやggplot2みたいに素人でも簡単便利だといいんだけど
Pythonはオブジェクト志向ぽいしすぐ諦めそう

ブラウザ環境なくなっちゃったの？

>>335
誤爆？
何をブラウズする環境が無くなったの？

ブラウザでプログラミングする環境

>>337
RStudio Serverの話？

RStudio Severあるで。Dockerで使うのがいいんじゃないかな？！

統計学とウェブ解析を交えて実践的な勉強と練習を
したいのですが、おすすめな書籍やサイトはありますか。
実際に解析ツールや分析ツールを用いて
自分で分析解析してから
解答を見て適切な手順や方法、考察を
解説してくれるものが良いです。


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

ウェブ解析とは具体的になんですか

R studioは日本語コメント書く度にIMEが無効になったりカーソルがずれたり黒文字の予測変換が黒背景と重なって見えなくなったりと散々だわ

Windows環境で使うから

会社規定なんでしかたない
UNIX環境使えるのが羨ましい

Docker for WindowsでRStudioサーバー動かせば？

なぜPythonスレッドがない？

板名読める？

>>342
IME無効になるのは俺だけじゃなかったと知ってほっとした。

ストアアプリも同じ症状でるからRStudio固有の問題でなくWindows環境の不治の病だと思ってる

Run any R code you like. There are over three thousand R packages preloaded.

https://rdrr.io/snippets/

>>350
日本語あると動作しないが、
グラフを描いてくれるのはうれしい。


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

Fisher test検定時に
p<2.2e-16
と表示されるんですが、これより小さい値の指数桁数を正確に表記する方法教えて下さい。
例えば5.8e-35となるようにです。

返値の中にあるp値を参照しなされ

fisher.test関数の返り値はlist型で、その中にp.valueという名前でp値が格納されているから$演算子を使って直接参照するか、broom::tidy関数に返り値を渡してdata.frame形式で出力してやれば見れる

>>352
>>354が言う直接的な参照
> fisher.test(matrix(c(1,120,130,2),2))$p.value
[1] 1.691912e-69

352です。よく分かりました！
ありがとうございます！

>>356
技術的な助言をしたけど、学術的に言えば、
p < 0.01 は全て p < 0.01 として、具体的なp値を考える意味はないと思うよ。
一部の例外的な研究分野を除いて(e.g., 遺伝統計学)。

はい、まさにその例外的な分野で使おうとしてます。ありがとうございます。

￥

￥

￥

￥

￥

￥

￥