【R言語】統計解析フリーソフトR 第6章【GNU R】 [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
2^127-1 =170141183460469231731687303715884105727は素数
Mathematicaは正しく計算してくれるけど
Rだと
> 2^127-1
[1] 170141183460469231738248242066868668888
と表示された。
あまり大きな数字は扱えないようだ。 2^1279-1も素数
10407932194664399081925240327364085538615262247266704805319112350403608059673360
29801223944173232418484242161395428100779138356624832346490813990660567732076292
41295093892203457731833496615835504729594205476898112116936771475484788669625013
84438260291732348885311160828538416585028255604666224831890918801847068222203140
521026698435488732958028878050869736186900714720710555703168729087 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 荒らし(◆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潰しどころか数学板さえ潰すことの出来ないでいる希ガスる。
作戦倒れしている希ガスる。
でも本人はなんとも思ってないんだろうな。哀れ。哀れ。哀れ。 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
参考資料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年ほど前にやってたな。
生命表と幾何ブラウン運動に基づく株価のモンテカルロシミュレーションで、今の貯蓄で死ぬまでにお金を使い切る確率を計算するの。 >>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 >>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()のヘルプを参照。 順番を表示させる、*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" >>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
お知恵を拝借したい。 自分なら、こんな感じ
> 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,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 ■ このスレッドは過去ログ倉庫に格納されています