探検
【R言語】統計解析フリーソフトＲ 第6章【GNU R】 [無断転載禁止]©2ch.net
レス数が1000を超えています。これ以上書き込みはできません。
?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)))
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))
}
実行すると
> 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を参照。
なぜわざわざ簡易生命表を使うの？
5～6年ほど前にやってたな。
生命表と幾何ブラウン運動に基づく株価のモンテカルロシミュレーションで、今の貯蓄で死ぬまでにお金を使い切る確率を計算するの。
公表されていればurl希望
第22回生命表（完全生命表）
ttp://www.mhlw.go.jp/toukei/saikin/hw/life/22th/index.html
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))
}
最新版の平成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例目が該当者であった。
この学年のスカトロ嗜好数の期待値を求めよ。
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)
2.5% 5% 50% 95% 97.5%
16 22 70 97 98
mean
66.31167
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)
coloroutパッケージ入れて、色を割り当てるとできるかも。
アーカイブに入っちゃっているけど、使える。
setOutputColors256()のヘルプを参照。
> n=21
> paste(n,ifelse(0<n%%10 && n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th'),sep='')
[1] "21st"
それだと、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"
なんでそこで && を使うかなあ。&& の意味わかってる？
ifelse には & を使うことになるのが普通だし、そうでなければ if ～ else のほうがよい。
ご指摘とパッケージの御教示ありがとうございました。
?&で勉強になりました。
御教示ありがとうございました。
toOrdinal のコードをみたら 11,12,13は別扱いしてあった。
if (tmp %in% c("11", "12", "13"))
tmp.suffix <- "th"
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)
同じ結果だよ
> 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)
結果は同じ
簡潔なコードありがとうございます。
AもBも賞金の期待値がどちらも10万なので
Aを何本買おうが合計の期待値は同じかなと思いつつ、sample関数を使ってシミュレーションしたのだけれど
一様分布にはならないのでシミュレーションの間違い、
もしくは自分が理解していないsample関数の特性を検出しただけなのか、と思っておりました。
rbinomでのシミュレーション結果とヒストグラムが一致して、一様分布になるという思い込みが間違いだと確信できました。
勉強になりました。ありがとうございました。
大したことではないけど
0から始まるので本数とindexを補正して
b <- function(){which.max(sapply(0:100, a))-1}
の方が正確だと思います。
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
と買い方に無関係にも思えるんだなぁ。
>合計の期待値の標準偏差は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
ご指摘ありがとうございます。
最大値を与える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)
とすると1000を超えているから同値の最大が複数あったことがわかりました。
ありがとうございました。
このスレは勉強になるなぁ。
hist(b2,breaks=20,col=sample(colors(),2))
とすると毎回配色の変わるツートンカラーのヒストグラムがでて面白い。
125です。レス遅くなったけど、やりたいことができました。ありがとう。
を，3行2列おきに平均して
,,1
[,1] [,2]
[1,] 5 17
[2,] 8 20
,,2
[,1] [,2]
[1,] 29 41
[2,] 32 44
としたいんですが，効率的な方法はないでしょうか。
実際のデータは，10000x10000x400なんです。
>3行2列
さっぱりわからん。
例えば、結果の[1,1,1]の5は、どれとどれの平均？
6行4列の左上の6つ1,2,3,7,8,9の平均5だと思う。
効率的かどうかはわからんが、やってみた
(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))
>>154です。
説明が悪くてすみません。
>>156の通りです。
ありがとうございます。
手元にデータがないので明日試してみます。
> 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
> 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
でよくね？
> 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
（変数名の重複を避けるため>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
>
上級者のコードを読むのは勉強になります。（深謝）
巨細にわたってご教示・ご議論いただきありがとうございました。
>>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氏のスキルに改めて感服。
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
> 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
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)
>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歳なので、老衰にせずに
主治医指定の病名で診断書を作成。
カルテ記載のしっかりした患者の診療は負担なくこなせて（・∀・）ｲｲ!!
剰余系での加減乗除での実際
> 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
そりゃあ、outer に渡す関数がまずいんだろ。
outer に渡す関数は、完全にベクトル対応になっていなければならないんじゃないの？
でなければ、わざわざ outer 使うこともなかろう。
レスありがとうございます。
自作関数から自作関数を呼び出しているのが原因かと思っていましたが、
ベクトル対応していないとはこういうことでしょうか？
加減乗は引数のどちらをベクトルにしても動作
> 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
outer から呼び出されたときの引数を見てみればわかるよ。たぶん1度しか呼ばれず、できあがりの行列の要素すべてを計算するための引数がベクトルで渡されているはずだから。
ありがとうございます。
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, ...)
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
ディープな世界やでほんま
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
$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
これってバグなのか？
> 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
ボロノイ図を制作した場合に分割された面積の平均を出力するには、どうすればよいのでしょうか？
宜しければどなたかお力添えをいただけないでしょうか？
SpatialPolygonsクラスまでできるといると仮定して、
各ポリゴンの中にareaスロットがあるから、
それを抽出して、その平均を取れば良い。
ボロノイ図がSpatialPolygonsクラスになっていないなら、
まずは変換。
描画領域の総面積÷点数 じゃないの？
何か暗黙の前提があるの？
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
と入力しボロノイ図は出せているのですが、
function (object, ...)
UseMethod("summary")
<bytecode: 0x06b5eef0>
<environment: namespace:base>
と表示されます。どの様に入力または再入力すればdir.areaが表示させられるでしょうか？
私自身初めてで、全く分からない状態でRを使用しているので、宜しければ分かり易く教えて頂けると幸いです。
検証用のデータがないと誰も追試しないと思う。
>154は模擬データがあったのでやる気になった。
例えばcor関数のmethodオプションの取り得る値は"pearson", "kendall", "spearman"みたいに。
そういうのは無いのかな？
初めて使う関数のオプションいちいちググるの面倒くさい。
ヘルプ表示じゃ駄目なのかい？
> args(cor)
function (x, y = NULL, use = "everything", method = c("pearson",
"kendall", "spearman"))
NULL
前にヘルプで見ようと思ったら書いてなかった覚えが…。マニアックな関数だったんで、そんなもんかな～と諦めてた。やっぱ基本はヘルプか。
>>277
argsは見れる関数と見れない関数があるのであんまり役に立たない印象。
>argsは見れる関数と見れない関数があるのであんまり役に立たない
だったら、例が悪いよ。見られない関数の例をあげないと
> ヘルプで見ようと思ったら書いてなかった
こっちだってわざわざ下請け関数のオプションまでは書かないよ。
ヘルプを見て未解決ならソースを見てくれ。
で R のメモリーをクリア
以前は「summary」としか入れていなかったのですが、「dd$summary」と入力したところ
全体のデータが出力されました。それ以後「dd$summary$dir.area」等で出せるようになりました。
大変初期的な場面で躓いておりましたが、現在何とか進んでおります。
色々とお答えくださった方々、とても参考になりました。本当にありがとうございました。
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
案外、高い確率になった。
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()
単なる二項分布。
コインが５回続けて表がでたら、0.5^5 <0.05なのでイカサマコインといわれちゃいそうなんだが、
１００回やってみると案外、５回表が続くので確率を計算しようと思ったが、解析的にできる頭がないので
シミュレーションしてみた。
> 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はここ
https://cran.r-project.org
プロジェクトの方はESETの誤検知っぽいんだよな
してねぇよｗｗ
totalvirusで調べると一件引っかかるし、Redirector検知だから、しょうがないね
走らせるたびに色がちがっておもしろい。
hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
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(1,-5,6)
でrunすると3じゃなくて6を返すんだけど、どこが間違ってる？
かっこが足りなかった
環境構築からもう面倒なイメージ
3系は数が全て小数扱いなので楽
Anacondaというパッケージでインストールすれば、今流行りのJupyter Notebookという開発環境で対話的にコーディングできる（Rも使える）
dplyrやggplot2みたいに素人でも簡単便利だといいんだけど
Pythonはオブジェクト志向ぽいしすぐ諦めそう
誤爆？
何をブラウズする環境が無くなったの？
RStudio Serverの話？
したいのですが、おすすめな書籍やサイトはありますか。
実際に解析ツールや分析ツールを用いて
自分で分析解析してから
解答を見て適切な手順や方法、考察を
解説してくれるものが良いです。
統計学は統計検定2級の知識はありますが
ウェブ解析はテキスト読んだだけです。
UNIX環境使えるのが羨ましい
IME無効になるのは俺だけじゃなかったと知ってほっとした。
https://rdrr.io/snippets/
日本語あると動作しないが、
グラフを描いてくれるのはうれしい。
http://imagizer.imageshack.com/img923/2879/m59LVv.png
p<2.2e-16
と表示されるんですが、これより小さい値の指数桁数を正確に表記する方法教えて下さい。
例えば5.8e-35となるようにです。
>>354が言う直接的な参照
> fisher.test(matrix(c(1,120,130,2),2))$p.value
[1] 1.691912e-69
ありがとうございます！
技術的な助言をしたけど、学術的に言えば、
p < 0.01 は全て p < 0.01 として、具体的なp値を考える意味はないと思うよ。
一部の例外的な研究分野を除いて(e.g., 遺伝統計学)。
多重比較だと意味あるかも
> 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
可能でしょうか？
例えばtemp <- c("0時間","8時間","24時間","48時間")として、
matplot()の引数にtempをとるやり方です。
他にもやり方あれば教えてください。
matplot(..., xaxt="n")
axis(1, at=seq(along=temp), lab=temp)
遅くなりましたがありがとうございました。
できました！
これどうなの？
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)
}
実態調査か何か？
<-と=は挙動が違う場合があるので、使い分けていますが、
代入はどっちかと問われたら、無論 <- または ->
なお、
> 1 -> x
これはエラーにならないけど、
> 1 = x
1 = x でエラー: 代入の左辺が不正 (do_set) です
これはエラー
俺は基本=派。
関数の定義は
z.test <- function(x,n=16,sigma=1){
z=sqrt(n)*mean(x)/sigma
2*pnorm(abs(z),lower=FALSE)
}
と書いている。
当たり前
そんな使い方なんてするかよ
他言語と同じく=一文字の方がすっきりしてイイ
x <- 1
if(x <- 1) print('YES')
if(x < -1) print('YES')
<-良くない
まず、= の場合
> z.test <- function(x,n=16,sigma=1){
+ z=sqrt(n)*mean(x)/sigma
+ 2*pnorm(abs(z),lower=FALSE)
+ }
> z.test(1:3)
[1] 1.244192e-15
<- で初期値設定すると、エラー
> z0.test <- function(x,n<-16,sigma<-1){
Error: unexpected assignment in "z0.test <- function(x,n<-"
> z=sqrt(n)*mean(x)/sigma
Error in mean(x) : object 'x' not found
> 2*pnorm(abs(z),lower=FALSE)
Error in pnorm(abs(z), lower = FALSE) : object 'z' not found
> }
= でも通常に動作する。
> z.test = function(x,n=16,sigma=1){
+ z=sqrt(n)*mean(x)/sigma
+ 2*pnorm(abs(z),lower=FALSE)
+ }
> z.test(1:3)
[1] 1.244192e-15
<- 推奨の人に聞きたいのだけど
<- でないと動作しないってことあるのだろうか？
=より面倒
うっとおしい
それと引数の指定は代入じゃないと思うんだが感覚が違うよかな？
essだとアンダースコアを入れるを勝手に「<-」になる。
>>394
それは代入ではなく、関数の規定値の設定。
規定値の設定は「=」と決められているので、「<-」はアウト。
ただし、関数を実行するときには「<-」を使うことができる。
> mean(x<-1:10)
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10
ありがとうございます。
そんな使い方ができたのですね。
こんなのしか知りませんでした。
> f <- function() {
+ x<<-1:10
+ mean (x)
+ }
> f()
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10
Hit <Return> to see next plot:
と出るの鬱陶しいな、と思っていたら
par(ask=FALSE)
と設定しておけばいいんだな。
data(netmeta::parkinson)
でなくて
data(parkinson, package = 'netmeta')
とするんだな。
だと、1色だけど
hist(rnorm(100),col=sample(colours(),(sample(1:10,1))))
なら1～10色で表示される。
何色に増やしてもいいがな。
何をやりたいのかな？
hist(rnorm(5000),col=apply(matrix(runif(80),4), 2, function(x){rgb(x[1],x[2],x[3],x[4])}))
ソース読みたくて
library(BayesFactor)
ttestBF_indepSample
では表示されなかったが、
library(BayesFactor)
BayesFactor:::ttestBF_indepSample
だと出てきた。
とりあえずこれ
https://github.com/sieste/STAN_cheatsheet
:::と3個なんだよなぁ
グダグダ言っていないで、さっさとヘルプを参照すればいいじゃないか？
?':::'
ありがとう。
stats:::t.test.default
でｔ検定のソースが見られた
例)
rm <- 1
としてしまった場合にrm関数を元に戻したいです。
よろしくお願いします。
直後ならRstudioでctrl+zで元に戻す。
再起動してたらこの方法は無理だろな
自分で作ったのではない関数の rm だったら rm(rm) で数値ベクトルオブジェクトが消えて、関数の rm が見えるようになる。
> rm <- 1
> ls()
[1] "rm"
> rm(rm)
> ls()
character(0)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]
深く考えなくてもrm(rm) で元に戻るけど。
rmが数値ではなくて自作関数だったら、少々ややこしいけど、
base::rm()で大丈夫だろう
既に説明されてるけど上書きじゃなくて別オブジェクト扱いになるからrmで消すか明示的にパッケージを::演算子で関数を指定してやればおけ
自作関数だと上書きされちゃうけど
> rm = function(x) asin(sqrt(x))
> rm(1)
[1] 1.570796
> rm
function(x) asin(sqrt(x))
> rm = 1
> rm
[1] 1
> rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
inherits = FALSE)
.Gloalenv（だったかな？）に自作関数があり、そこに上書きして変数にしてるから消したら消えるだけで戻せない
（Rの）環境問題について その１。
ttps://qiita.com/kohske/items/325bdf48f4f4885a86f1
戻せます。
> rm <- function(x) asin(sqrt(x))
> rm
function(x) asin(sqrt(x))
> base::rm(rm)
> rm
function (..., list = character(), pos = -1, envir = as.environment(pos),
[以下略]
>.Gloalenv
case sensitiveなRでそれはないのでは？
> help(".GlobalEnv")
みなさまご返信ありがとうございました。
例示したものがbase::rmだったので自作関数のことまでは想像だにしませんでした…(がとても良い勉強になりました。)
自環境でも上書きした関数を元に戻すことが出来ることを確認しました。
戻したいのが自作関数
rm = function(x) asin(sqrt(x))
だったら上書きされて無理じゃないの？
ゴミ箱から消去したファイルを復活させるソフトもあるから
PCに詳しければ可能かもしれないけれど
Rの機能としては復活させるのは無理と思う。
base::rm(rm)で自作の方のrmは消えてbase::rmが「復活」します。
自作の方を上書きした場合は、そっちの復活は無理ですね。
ちょっと意味がわからない。
自作の関数なら定義を再実行するだけなのに、
消える消えないって何が問題なんだ？
関数だけでなく既定値でも同じ。
> pi
[1] 3.141593
> pi=2
> pi
[1] 2
> rm(pi)
> pi
[1] 3.14159
https://seleck.cc/685
https://trello.com/
ボードのメニュー → Power-Upsから拡張可能 Slack DropBoxなど
Trello Chrome拡張機能 elegant
ttp://www.kikakulabo.com/service-eft/
trelloのオープンソースあり
共同ツール 2
https://www.google.com/intl/ja_jp/sheets/about/
共同ツール 3
https://slack.com/intl/ja-jp
https://www.dropbox.com/ja/
https://bitbucket.org/
https://ja.atlassian.com/software/sourcetree
https://www.sketchapp.com/
ttp://photoshopvip.net/103903
ttps://goodpatch.com/blog/sketch-plugins/
Trello Chrome拡張機能プラグイン集
https://chrome.google.com/webstore/search/trello?_category=extensions
Slackプラグイン集
https://slack.com/apps
Sketchプラグイン集
https://sketchapp.com/extensions/plugins/
https://supernova.studio/
知恵袋に漸化式があったので計算スクリプトを書いてみた。
表がでる確率ｐのコインをN回投げてK回以上表が続く確率。
# N: flips
# K: least sequential head
# p: probability of head
seqNp <- function(N=100,K=5,p=0.5){
if(N==K) return(p^K)
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q for any i
for(i in K:(N-1)){ # recursive formula
a[i+1]=0 # a[i+1]=q/p*(a[i]+a[i-1]+a[i-2]+...+a[i-(K-1)])
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}
P0=numeric(N) # P0[n] : probability of ending with 0 head when flipped n times
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0
# MP[k,n] = Pk[n] : probability of ending with k head when flipped n times
for(k in 2:K){ # Pk(n+1)=p*P(k-1)(n)
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
}
ret=1-apply(MP,1,sum)
ret[N]
}
この例で、答えとして何が返ってきて欲しいのか？
2で合ってる？
1が5回以上出現しているからTRUE（数字なら１）が返り値。
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.
chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}
whichとdiffを用いての解法のお返事ありがとうございます。
実行時間を>302のスクリプトと比べてみました。
run.check <- function (x, n=5) {
# x の中に 1 が n 回以上連続していれば TRUE を返す.
chg <- c(TRUE, diff(x) != 0) # 変化があった場所
chgidx <- c(which(chg), length(x)+1) # 変化があった場所の添え字
run.length <- diff(chgidx) # 0や1の連続している個数
true.length <- run.length[x[chg] == 1] # 1の連続している個数
any(true.length >= n) # 連続している個数が n 以上のrunがあるか?
}
# N(=100)回コインをなげてｎ(=5回)以上続けて表がでる確率。
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
> system.time(mean(replicate(10^5,run.check(rbinom(100,1,0.5)))))
user system elapsed
17.56 0.04 17.68
> system.time(mean(replicate(10^5,seqn())))
user system elapsed
9.74 0.07 9.78
後者の方が速いのは1がn個連続したら、TRUEを返して終了するので以後のチェックはしないためだろうと思います。
diff や anyの使い方が勉強になりました。
時間を割いてスクリプトを作成していただいてありがとうございます。<(_ _)>
具体的には、IME自体をオンにすることができません。
日本語の表示やコピペは問題なくできます。
環境
Linux - Utubntu
RStudio Version 1.0.143
パッチはUbuntu 16.04.LTS用にしかないので他のバージョンならRStudioサーバー使った方が早いかと
なお、パッチはRStudio 日本語とかでグクってみて
別のソフトでIMEのON/OFFをやってから戻ると直るけど。
うちも同じ。
Ubuntuはfcitxのパッチ当てとけばWindowsのような症状は出ないので快適
http://img-cdn.jg.jugem.jp/d75/1405116/20140609_1105745.jpg
スクリプトは
http://imagizer.imageshack.com/img921/8020/t1VcqA.jpg
小数点以下500桁でみると19個素数があった。
> vip=Vectorize(is.prime)
> (idx=which (vip(N)))
[1] 100 124 150 172 183 202 215 219 255 296 310 314 323 346 400 417 439 441 474
> cat(substring (e,idx[1],idx[1]+9))
7427466391
９９人の囚人がいます。彼らの頭に１～１００までのナンバーカードが貼りつけられた帽子をランダムにかぶせます。
他人の帽子は見ることができても、自分の帽子は見ることができません。
帽子の数は全部で１００なので、一つ使われずに余ります。
そのナンバーは囚人達にはわからないようにしておきます。
この状況で、囚人たちに一斉に自分のナンバーを宣言させて、全員が正解だったら釈放するという賭けをします。
囚人たちには帽子をかぶせられる前に相談タイムが設けられています。
どういう戦略を取れば、助かる確率を最も高くできるでしょうか？
http://000013.blogspot.com/2010/12/99.html
解答を読んでも数理が理解できないが
確率が0.5になるのはシミュレーションできた。
# http://000013.blogspot.com/2010/12/99.html
inversion <- function(x){
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret) # inversion number
}
is.even= function(x) !inversion(x)%%2 # is inverion number even?
prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx])
Y=numeric(n-1)
for (i in 1:(n-1)){ # select as even permutation
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # two numbers unseen for i-th prisoner
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2])
}
all(X==Y)
}
mean(replicate(1e3,prisoner99()))
遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ
どうもＣでやったらうまくいかないのでＲでやってみた。
# RED + WHITE = COLOR
x=c('R','E','D','W','H','I','T','E','C','O','L','O','R')
unique(x)
redwhite <- function(R,E,D,W,H,I,T,C,O,L){
red=10^(2:0)
white=10^(4:0)
color=10^(4:0)
sum(red*c(R,E,D))+sum(white*c(W,H,I,T,E)) - sum(color*c(C,O,L,O,R))
}
x=unique(x) ; x
REDWHITECOLOR <- function(x){
R=x[1]
E=x[2]
D=x[3]
W=x[4]
H=x[5]
I=x[6]
T=x[7]
C=x[8]
O=x[9]
L=x[10]
cat(paste('RED = ',R,E,D,
' WHITE = ',W,H,I,T,E,
' COLOR = ',C,O,L,O,R),'\n')
}
REDWHITECOLOR(unique(x))
library(gtools)
perm=permutations(n=10,r=10,v=0:9)
perm=perm[perm[,1]!=0&perm[,4]!=0&perm[,8]!=0,] # R!=0,W!=0,C!=0
head(perm) ; tail(perm)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=redwhite(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6],perm[i,7],perm[i,8],perm[i,9],perm[i,10])
}
hist(re)
(indx=which(re==0))
REDWHITECOLOR(perm[indx,])
> REDWHITECOLOR(perm[indx,])
RED = 5 8 7 WHITE = 3 9 6 1 8 COLOR = 4 0 2 0 5
# × CD
# ------
# CCE
# BFAA
# ------
# BEBDE
f <- function(A,B,C,D,E,F){
(A*100+A*10+B)*D==C*100+C*10+E &(A*100+A*10+B)*C==B*1000+F*100+A*10+A &
(C*100+C*10+E)+(B*1000+F*100+A*10+A)*10==B*10000+E*1000+B*100+D*10+E
}
library(gtools)
perm=permutations(n=10,r=6,v=0:9)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=f(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6])
}
indx=which(re==TRUE)
perm[indx,]
PERT分布を学んでいるのですが、ベータ分布のスケーリングの仕方が分かりません。
最小値a、最頻値ｂ、最大値ｃ として、a, b, c から求めたパラメータをa1、a2 とします。
このとき、
dbeta(x, a1, a2)*(c-a)+a
とすると、当然ながら横軸の値の範囲は0から1で、縦軸の値のみ大きくなってしまいます。
横軸の値の範囲がaからcで、縦軸の値は確率密度のままにするには、どうしたらよいのでしょうか？
dbeta(x-a, a1, a2)
447です。ご教示どうもありがとうございました。
a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
mu <- (a+4*b+c)/6 # PERT分布の平均値
a1 <- b*(mu-a)/(c-a) # パラメータ a1
a2 <- b*(c-mu)/(c-a) # パラメータ a2
pert1 <- function(x) dbeta(x,a1,a2) # ベータ分布布の計算
まではできたのですが、ご教示いただいた方法で
pert1 <- function(x) dbeta(x-a ,a1,a2)
とすると、ｘ の範囲が0から1で、ｙが0の一様分布のようなグラフになってしまいます。
https://en.wikipedia.org/wiki/Beta_distribution#Four_parameters_2
a1とa2のパラメータ推定って合ってるんだっけ？僕も正直わかんない
やりたいことはこれだと思う
install.packages("mc2d")してからやってね
----------------------------------------------------
rm(list=ls())
library(mc2d)
a <- 40 # 最小値
b <- 45 # 最頻値
c <- 70 # 最大値
pert4<-function(x) dpert(x,min=a,mode=b,max=c,shape=4)
x1<-40:70
plot(x1,pert4(x1))
447です。
450様、まさにこれ、これです。どうもありがとうございました。
おかげさまで大変助かりました。
448様にも、ご助言下さいましたこｔ、厚く御礼申し上げます。
Windows 10 （Rの日本語ヘルプはインストールしていない）+ R studioだけど、そのメッセージ出ないな。
「安かろう悪かろう」「俺が知らないものは三流」的な考え方の人が多くて肩身が狭い…
>>452のようなエラーが出たら「無料のものはやっぱりダメね」
有償(SPSS) →信頼できる
ボランティア→いい加減
と言う思考らしい
ユーザー名はアルファベット1文字です
studioインストールしてないんですけど、しないとダメですか？
> Sys.setlocale(locale="C") # locale の変更（USの設定にしないと表示が変になる）
[1] "C"
元に戻すには
> Sys.setlocale(locale="Japanese_Japan.932")
フィッシャーの正確確率検定で、
適合度（一様性）の検定ってできますか？
できるのは独立性の検定のみ？
Rでいろいろ試したみたんですが、
2群ないとエラーになる。。
レス感謝です。
そうなんですが、
期待度数が5未満になってしまった場合は、
カイ二乗検定は向かないんですよね？
こういうので
r1=5;r2=4;n1=10;n2=12
prop.test(c(r1,r2),c(n1,n2),correct=TRUE)
chisq.test(matrix(c(r1,r2,n1-r1,n2-r2),2,byrow=TRUE),correct=TRUE)
警告がでるという話かな？
Warning message:
In chisq.test(matrix(c(r1, r2, n1 - r1, n2 - r2), 2, byrow = TRUE), :
Chi-squared approximation may be incorrect
適合度を見るのに関係ないだろ。
> cnt
[1] 1 0 2 4 3 4 12 6 12 11 21 22 21 37 40 30 59 44 49 47 38 48 36 43 38
[26] 46 33 34 30 21 26 33 17 13 14 12 8 16 7 7 9 6 6 5 3 8 4 4 4 1
[51] 0 2 1 0 0 0 0 2
> prb
[1] 0.0002947255 0.0006464052 0.0012538873 0.0022037831 0.0035720843
[6] 0.0054114364 0.0077413658 0.0105430123 0.0137588816 0.0172972083
[11] 0.0210398698 0.0248524558 0.0285950600 0.0321325300 0.0353432214
[16] 0.0381256526 0.0404028045 0.0421240952 0.0432652773 0.0438266419
[21] 0.0438299769 0.0433147339 0.0423338191 0.0409493611 0.0392287274
[26] 0.0372409836 0.0350539129 0.0327316483 0.0303329196 0.0279098757
[31] 0.0255074215 0.0231629878 0.0209066521 0.0187615256 0.0167443293
[36] 0.0148660914 0.0131329065 0.0115467110 0.0101060386 0.0088067297
[41] 0.0076425765 0.0066058952 0.0056880193 0.0048797146 0.0041715201
[46] 0.0035540195 0.0030180499 0.0025548579 0.0021562071 0.0018144483
[51] 0.0015225575 0.0012741483 0.0010634658 0.0008853653 0.0007352804
[56] 0.0006091856 0.0005035529 0.0004153084
> chisq.test(cnt, p=prb, rescale.p=TRUE, simulate.p.value=TRUE)
Chi-squared test for given probabilities with simulated p-value (based
on 2000 replicates)
data: cnt
X-squared = 65.504, df = NA, p-value = 0.2189
5未満の数があっても関係がない
警告が出るのは独立性検定の場合でしたね。
カイ二乗検定で、期待度が5未満が多いのは不適切なのかと思い込んでおりましたが、
それは独立性の検定の場合のことであり、
適合度（一様性）の検定の場合は、
気にしなくて良かったのか。。
ならカイ二乗検定でやります。
いろいろ教えて下さり、ありがとうございます。
http://www.cdleidyba.lt/file/spss_old_1.png
Stanとかに対応してんの？
化石みたいなバージョンやな
いまSPSSを使う理由って何なの？煽りじゃなくマジで
周囲のspssユーザを見ると「指導教員がspssしか使えないので、
自分もspssしか使えない、
新しい（学術系）ソフトは誰かに丁寧に根気よく指導してみらえないと無理」、
これの再生産。
研究室に予算がないから
予算がなかったら、Rを使え
フィボナッチ数列の再帰呼び出し
# 関数名での呼び出し
fibo <- function(n){
if(n==1|n==2) return(1)
else fibo(n-1)+fibo(n-2)
}
# Recallを使っての呼び出し
fiboR <- function(n){
if(n==1|n==2) return(1)
else Recall(n-1)+Recall(n-2)
}
> system.time(fibo(30))
user system elapsed
4.27 0.02 4.36
> system.time(fiboR(30))
user system elapsed
6.65 0.00 6.70
何度も試行錯誤
# cumsum with recursive call
cumsumR <- function(v){
res=numeric(length(v))
cumsumR_sub <- function(v,res,i){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
cumsumR_sub(v,res,1)
}
cumsumR(1:10)
とりえあず、動作した。
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
cumsumR <- function(v,res=NULL,i=1){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
>
これってifelseの仕様でしょうか？バグでしょうか？
> x=2:1
> if(x[1]<=x[2]) x else x[2:1]
[1] 1 2
> ifelse(x[1]<=x[2],x,x[2:1])
[1] 1
ifelse() はベクトル演算できる関数。
だから、戻り値の要素数は、test部分の要素数に一致する。
testが1つなので、x[2:1]の1つ目の要素が返される。
> ifelse(c(x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2
さらに３つにすると、helpに書いている通り、yesやno部分は再利用される。
> ifelse(c(x[1] <= x[2], x[1] <= x[2], x[1] <= x[2]), x, x[2:1])
[1] 1 2 1
ありがとうございました。
バグではなくてそういう仕様だったのですね。
こういう使い方ができるということがわかって勉強になりました。
> x=2:1
> ifelse(c(x[1]!=x[2],x[1]==x[2]),1:2,3:4)
[1] 1 4
ifelse() はベクトル演算できると、再利用されるの動作確認。
> ifelse(c(TRUE,TRUE,FALSE,FALSE,TRUE,FALSE,FALSE),1:2,11:15)
[1] 1 2 13 14 1 11 12
> mat
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 0
[3,] 1 1 0
[4,] 1 1 0
[5,] 1 1 0
[6,] 1 1 0
> print(apply(mat,2, function(x) ifelse(x==1,'真','偽')),quote=F)
[,1] [,2] [,3]
[1,] 真 真 真
[2,] 真 真 偽
[3,] 真 真 偽
[4,] 真 真 偽
[5,] 真 真 偽
[6,] 真 真 偽
> (mat=matrix(sample(0:1,12,rep=T),3,4))
[,1] [,2] [,3] [,4]
[1,] 1 1 0 1
[2,] 0 1 0 0
[3,] 1 0 0 1
> print(ifelse(mat==1,'Right','Wrong'),quote=F)
[,1] [,2] [,3] [,4]
[1,] Right Right Wrong Right
[2,] Wrong Right Wrong Wrong
[3,] Right Wrong Wrong Right
>
エラー処置は省略。
reduce_fraction <- function(x,y){
a=x
b=y
r = a%%b # a=bq+r -> a%%b=b%%r
while(r){
a = b
b = r
r = a%%b
}
gcd=b
cat(x/gcd,'/',y/gcd,'\n')
invisible(c(x/gcd,y/gcd))
}
> reduce_fraction(2860,1082900)
11 / 4165
login_page <- html_session("https://xxxxx";)
login_form <- html_form(login_page)[[1]] %>%
set_values(AAA="xxxx", BBB="xxxx")
というので行けるはずなんですが
最終行のAAAとBBBってそれぞれhtmlタグのname=""のとこからとってくるんですよね?
これにここがname="user[email]"みたいに[]という記号はいってたらどうすればいいでしょうか?
よろしくおねがいいしまーーーーーーーす
質問の意味がいまいち分からないから、誰も助言できないのでは？
> x <- 'user[email]'
> x
[1] "user[email]"
rvestパッケージを使ったことがないけど、
記号が入っていたらなぜ問題があるのかよく分からない。
仮でもいいから
ちがいます
AAAがどれかよくみてみてください
””はありません
<input class="form-control" autofocus="autofocus" placeholder="Email address" required="required" type="text" value="" name="user[email]" id="user_email">>>494
もっと分かるように説明しないと、追試できる情報も提供しないし。
もしかして変数名などに[]が入っている場合 にどうしたらよいかってこと？
> `user[email]` <- 1:5
> `user[email]`
[1] 1 2 3 4 5
多分行けました
超感謝!!
問題文の記号のまま
q <- function(x)
....
とやって
q(100)と入力すると
Rが終了することに気づいた。
を実行すると爆弾出るね
空ベクトルの検出って長さ０以外に検出方法ってあるでしょうか？
> x=c(1,2)
> x=x[-1]
> x=x[-1]
> length (x)
[1] 0
> x
numeric(0)
> x==numeric (0)
logical(0)
> is.null(x)
[1] FALSE
> is.na(x)
logical(0)
> x==NULL
logical(0)
> x==NA
logical(0)
> length(x)==0
[1] TRUE
高さはパッケージnleqslvを使った近似計算。
# Calculate tetrahedron volume from cordinates
library(nleqslv)
Tetra <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
fn <- function(x,O,A,B,C){
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO # H on triangle ABC
AB=B-A
AC=C-A
c(HO%*%AB,HO%*%AC) # HO vertial to AB and AC
}
fn1 <- function(x) fn(x,O,A,B,C)
x=nleqslv::nleqslv(c(1/3,1/3),fn1)$'x'
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO
h=sqrt(sum(HO^2))
a=sqrt(sum((B-C)^2))
b=sqrt(sum((C-A)^2))
c=sqrt(sum((A-B)^2))
s=(a+b+c)/2
base=sqrt(s*(s-a)*(s-b)*(s-c))
V=1/3*base*h
return(V)
}
初期値は辺の長さ１の正四面体
> options(digits = 16)
> Tetra()
[1] 0.1178511301977579
> sqrt(2)/12
[1] 0.1178511301977579
多分、正常に動作していると思う。
行列式det使うと簡単
> po <- c(1/2, sqrt(3)/6, sqrt(2/3))
> pa <- c(0,0,0)
> pb <- c(1,0,0)
> pc <- c(cos(pi/3), sin(pi/3), 0)
> det(cbind(pa-po,pb-po,pc-po))/6
[1] -0.1178511
tetrahedron <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
AO=A-O
BO=B-O
CO=C-O
as.numeric(abs(pracma::cross(AO,BO) %*% CO)/6)
}
４行で済んだ。
>>503
ありがとうございました。
パッケージに頼らずに計算できたのですね。
私が思いつくのは
summary
boxplot
hist
pairs
です。こんなのも良いよってのがあったら教えてくださいm(_ _)m
※ライブラリの使用有無は問いません
BESTパッケージのplotPost
histに加えて95%CIを表示してくれる。
出力をデータフレーム型にしたいときはskimrパッケージ
良い方法が見つからなかったので !length(x)で空白ベクトル判定とした。
文字列を逆に並べる再帰呼び出しスクリプト
> reverse <- function(x){
+ if(!length(x)) return(NULL)
+ c(Recall(x[-1]),x[1])
+ }
> cat(reverse(LETTERS[1:26]))
Z Y X W V U T S R Q P O N M L K J I H G F E D C B A
なぜ if (length(x) == 0) と書かないのか？
可読性って知ってるか？
文脈読めば
!length(x)で空白ベクトル
の流れ
length は数値を返す関数なのに、なぜわざわざ論理演算をするんだ？
プログラムは、ここの議論を知らない人が読む可能性のほうが高いんだぞ。
正規分布が一様分布より大きい期待値の算出に
mean(rnorm(100)>runif(100))
って書かない？
俺は
sum(ifelse(rnorm(100)>runif(100),1,0))/100
って書いたりしないけど。
え、常套手段やん？！
論理値を数値に置き換えて計算しているから
数値を論理値にしても別に違和感がない。
可読性は慣れの問題。
アンカー間違えとる？
だから、そんなの常套手段やん？
それとこれとは話が別だ。
logical は TRUE か FALSE の二値しかとらないから、それから数値への変換は自明。
多種の値をとる数値に論理演算をするのはいただけない。
is.empty.vector のような関数があるならそれでよいが、length を使っているのだから、それを勝手にベクトルが空かどうかの判断に使っているのはあなたであって、他人はそのようには考えない、といっているのだ。
自分だけしかこーどを見ないならべつに構わないが、こういうところに晒すのはよくない。
ましてや常套手段などと言って他人に教え込むのはやめてもらいたい。
常套手段というのはmeanの話じゃね？
n=0
while(1){
if(n>10) return(10)
n=n+1
}
n=0
while(2){
if(n>10) return(10)
n=n+1
}
配列を逆順に並べる再帰呼出しのコード。
Rにはrevという関数がある。
そのコードみてみ！
lengthを真偽判定に使ってますがな。
> base:::rev.default
function (x)
if (length(x)) x[length(x):1L] else x
>521みた？
if(length(x)!=0)になってないよ。
>こういうところ
誰も鵜呑みにしない便所の落書き。。。
if(!0) print(!1)の結果もサクッとわかる。
関数によるな。
if やwhileは変換処理しているけど
whichだとエラーになるな。
> which(!0)
[1] 1
> which(0)
Error in which(0) : argument to 'which' is not logical
> which(2)
Error in which(2) : argument to 'which' is not logical
> which(!3)
integer(0)
>
> which(!(-1:1))
[1] 2
> !(-2:2)
[1] FALSE FALSE TRUE FALSE FALSE
is.zero だと？
だったら最初から == 0 てよいじゃないか。
だいたい、if や while は関数じゃないし。
while (1) などと書くのはCに毒されているんじゃねーの？
Rなら while (TRUE) のほうが良いし、repeat というのもある。
スクリプトを読むときに脳内変換する癖をつけておけば
可読性が悪いと思わずにすむ。
>521みた？
if(length(x)!=0)になってないよ。
Rのライブラリ？の書き方なんてあてにならないよ。
そもそも、rev も 509 の reverse も演算の回数を必要最低限にするなら、条件は length(x)) >= 2 または > 1 だ。
まあ if (length(x)) のほうが速かったのかも知れないが、常にそうなるとは限らないのだから、早すぎる最適化は諸悪の根源だ。
自分の流儀と違う { の使い方でも文句いいそう。
!
!=
>
==
での真偽判定を各々１０回する時間を出してみた。
> x=rbinom(1e7,1,0.5)
> system.time(replicate(10,!x))
user system elapsed
1.27 0.33 1.59
> system.time(replicate(10,x!=0))
user system elapsed
2.39 0.36 2.75
> system.time(replicate(10,x>0))
user system elapsed
2.58 0.31 2.92
> system.time(replicate(10,x==1))
user system elapsed
2.60 0.36 2.99
> f <- function(){
+ x=rbinom(1e7,1,0.5)
+ (a=system.time(!x)[3])
+ (b=system.time(x!=0)[3])
+ a<b
+ }
> (re=replicate(10,f()))
elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed elapsed
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
elapsed
TRUE
結果は再現された。
> x=1:1e8
> f1 = function(x) if(!length(x)) return(NULL) else x=x[-1]
> f2 = function(x) if(length(x)==0) return(NULL) else x=x[-1]
> system.time(f1(x))
user system elapsed
2.38 0.53 2.93
> system.time(f2(x))
user system elapsed
2.05 0.61 2.69
御指摘の通り、!の負け
＞遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
＞詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ
Ubuntu 18.04.1 LTSにアップグレードしたけどできない
いままでスクレイピングにつかってrvestのログインコマンドもエラー出るようになったし最悪だ。。。。
Voyager 18.04LTS（xubuntu 18.04LTS）ではパッチなしでRStudioに日本語入力できた。
作業ディレクトリの固定化できなくね？
前まではショートカットのプロパティで作業フォルダ指定してたんだけど、
新しいR、何度やっても作業ディレクトリがデフォルトのドキュメントに戻りやがる
どうすればいいか、誰がご教示を。
プロパティでコマンドラインオプションを消せばよい。
--cd-ほにゃらら
というのが付いてるでしょ。
あ、なるほど
ありがとうございます、多謝
以下のコードで、dataには2018/10/20 00:00～2018/10/25 00:00の1時間刻みのaとbの値のxts型オブジェクトが入りますが、（aとbは時系列データのイメージです）
このdataの毎日09:00の値を添え字を使って取り出す方法はありますか？
（data["2018-10-21 00:00"]とすると10/21の00:00が取り出せるような感じで）
library(xts)
time <- seq(as.POSIXct("2018-10-20 00:00"), as.POSIXct("2018-10-25 00:00"), by="1 hour")
a <- rnorm(length(time))
b <- rnorm(length(time))
data <- xts(x = cbind(a,b), order.by = time)
思いついたのはfor文を使って以下のようにするのかなと思ったのですがもう少しいい方法があるだろうなと思いまして…
data_result <- data[10] # まずxts型の入れ物を作る(2018/10/20 09:00)
for (i in 1:4) {
data_result <- rbind(temp,data[10+24*i]) # 2日目以降はこれにくっつける
}
data_result
data_result <- data[10] # まずxts型の入れ物を作る(2018/10/20 09:00)
for (i in 1:4) {
data_result <- rbind(data_result, ,data[10+24*i]) # 2日目以降はこれにくっつける
}
data_result
idx <- seq(10, length(time), by=24)
result <- data[idx]
じゃダメなのか？
ああ、これがいいですね！お恥ずかしい…
ありがとうございます
５０の階乗を計算させてみた
R with Rmpfr
> mpfr(factorial(50),1e5)
1 'mpfr' number of precision 100000 bits
30414093201713018969967457666435945132957882063457991132016803840
Haskell:
Prelude> product[1..50]
30414093201713378043612608166064768844377641568960512000000000000
Python:
import math
print(math.factorial(50))
30414093201713378043612608166064768844377641568960512000000000000
Wolfram:
https://www.wolframalpha.com/input/?i=50!
30414093201713378043612608166064768844377641568960512000000000000
enc2native()しても駄目だし
とりあえずcairo_png()だけでも…
ある医院に１時間あたり平均５人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
1時間あたりの平均診療人数は６人で、一人あたりの診療時間は指数分布に従うとする。
診察までの平均の待ち時間は何時間か？
このシミュレーション解はこれであってる？
N=1e6
λ=5
μ=6
sum(rpois(N,λ)*rexp(N,μ))/N
> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.8331036
60かけて、分にすると
[1] 49.86565
grep 使って書いてみた。
# mhs(c(1,0,1,1,0,1,1,1)) : return 3
mhs = function(x){ # maximum head sequence
y=paste(x,collapse='')
str="1"
if(!grepl(str,y)) return(0)
else{
while(grepl(str,y)){
str=paste0(str,"1")
}
return(nchar(str)-1)
}
}
system.time(mean(replicate(1e4,mhs(sample(0:1,100,rep=T))>=5)))
# N(=100)回コインをｎ(=5回)以上続けて表がでるか？TRUE or FALSE
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
system.time(mean(replicate(1e4,seqn())))
結果は、
> system.time(mean(replicate(1e4,mhs(sample(0:1,100,rep=T))>=5)))
user system elapsed
4.56 0.01 4.61
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
1.05 0.00 1.07
逆に４倍時間がかかるようになった
そのベクトルの5個前とのdiffの要素の中に5が一回でも現れることが、一回でも5回連続で表が出たことがあることの必要十分条件
grepはどうしても時間がかかるしif文もできれば使いたくない
searchseq <- function(n=5,N=100,p=0.5,trial=10000){
result <- 0 # 表が5回以上出た回数を数える入れ物
for (i in 1:trial){
coin <- rbinom(N,1,p)
coincumsum <- cumsum(coin)
coindiff <- diff(coincumsum,5)
#diff(coincumsum,5)の要素に一個でも5があれば表が5回以上出たことがあったということ
#anyでT/Fにして、sumで0/1にする
result <- result+sum(any(coindiff==5))
}
return(result/trial)
}
結果もよさそう
> searchseq()
[1] 0.80793
> system.time(mean(replicate(10000,mhs(sample(0:1,100,rep=T))>=5)))
ユーザ システム 経過
0.92 0.00 0.93
> system.time(mean(replicate(10000,seqn())))
ユーザ システム 経過
0.29 0.00 0.30
> system.time(searchseq())
ユーザ システム 経過
0.21 0.00 0.20
レスありがとうございました。
rleとも比べてみました。
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
4.290 0.000 4.377
>system.time(mean(replicate(1e4,max(rle(rbinom(100,1,0.5))$len)>=5)))
user system elapsed
3.620 0.000 3.742
> system.time(mean(replicate(1e4,mhs(rbinom(100,1,0.5)>=5))))
user system elapsed
2.390 0.000 2.435
> system.time(searchseq())
user system elapsed
1.880 0.000 1.988
greplの逆転はsampleをrbinomに変えたためでしょう。
> system.time(replicate(1e5,sample(0:1,100,rep=T)))
user system elapsed
7.200 0.180 7.591
> system.time(replicate(1e5,rbinom(100,1,0.5)))
user system elapsed
5.980 0.230 6.319
rleは0も数えているから間違っているんじゃ？
ご指摘ありがとうございました。
修正しました。
rle1 <- function(N=100,n=5,p=0.5){
r=rle(rbinom(N,1,p))
max(r$len[which(r$val==1)])
}
結果
> system.time(mean(replicate(1e4,max(rle1()>=5))))
user system elapsed
4.430 0.000 4.546
> system.time(mean(replicate(1e4,seqn())))
user system elapsed
4.490 0.010 4.609
> system.time(mean(replicate(1e4,mhs(rbinom(100,1,0.5))>=5)))
user system elapsed
7.440 0.000 7.656
> system.time(searchseq())
user system elapsed
1.950 0.000 2.066
>
system.time(mean(replicate(1e4,any(diff(cumsum(rbinom(100,1,0.5)),5)==5))))
user system elapsed
1.820 0.000 1.886
>
> system.time(mean(replicate(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[wh
<e(1e4,with(rle(rbinom(100,1,0.5)), max(lengths[whi ch(values==1)])>=5))))
user system elapsed
4.370 0.010 4.478
意味もなくforループ回してた上に毎回sum使って真偽値を数値に変換してたけど
replicate使って最後に一回だけmean取ると2.066→1.886で1割短くなるのね
他人のコード読むのは勉強になる
whichいらねーよ。
全部が０のとき、エラーになるので修正
rle01 <- function(x){ # c(0,1,1,1,0,0) => return 3
if(sum(x)==0) return(0) #c(0,0,0,0,0,0) => return 0
else{
r=rle(x) # Run Length Encoding
max(r$lengths[which(r$values==1)]) # max length of value 1
}
}
動作確認
> rle01(x<-rbinom(100,1,0.5)) ; x
[1] 8
[1] 1 1 0 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1 0 0 0 1 1 1 1 0 0 1 1 1 1 1 0 0
[36] 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 1 1 0 1 1 0 0 1 1 0 1 0 1 1 0 1 0 1 0 0
[71] 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 1 1 1 1 0 0 1 0 0 1 1 1 0 1 1
> rle01(x<-rbinom(100,1,0)) ; x
[1] 0
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[36] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[71] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
センスないなあ。
rle01 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}
whichは余分でした。
sumの方がrleより高速だと思ったらから、すべて０の場合はrleを呼ばないことにしただけ。
0連続でやるとこれだけ差がつく。
rle01 <- function(x){ # c(0,1,1,1,0,0) => return 3
if(sum(x)==0) return(0) #c(0,0,0,0,0,0) => return 0
else{
r=rle(x) # Run Length Encoding
max(r$lengths[which(r$values==1])] # max length of value 1
}
}
rle012 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}
> x=rep(0,1e8)
> system.time(rle01(x))
user system elapsed
0.3 0.0 0.3
> system.time(rle012(x))
user system elapsed
7.36 4.52 13.25
>
そんな特別な場合のことに対して高速化するのは愚の骨頂。
1が一つでもある場合にsumを呼ぶのは余計じゃないのか？
余計なwhichがあるくらいだから、まずは素直にやりたいこと、やるべきことを正しく書くようにしたら？
確率の理論値は9.8897353347449091e-05
> system.time(mean(replicate(1e4,rle01(rbinom(N,1,p))>=n)))
user system elapsed
0.61 0.00 0.64
> system.time(mean(replicate(1e4,rle12(rbinom(N,1,p))>=n)))
user system elapsed
1.97 0.00 2.03
分数で表すと
890788167367/9007199254740992
9.889735334744909e-05
どんな引数をつければいいですか
実はシミュレーションじゃなくて
漸化式からのプログラム解を分数表示するプログラムはpythonで作成済。
ここに置いた。
https://egg.2ch.net/test/read.cgi/hosp/1540905566/77
このコインの表がでる確率の期待値、モード比、信頼区間を求めるのが次のネタ。
unirootで算出できたけど
シミュレーションはどうすればいいのかアイデアが浮かばない。
MCMCで解決できるかなぁ？
これをシミュレーションで検証したい。
$Rscript main.r
lower mean mode upper
0.2487456 0.4469764 0.4589692 0.6386493
http://tpcg.io/asKRE9
# maximal sequential head probability at 10^8 coin flip
> y
[1] 2.2204460492503131e-16 2.2204460492503131e-16
[3] 8.8817841970012523e-16 1.5543122344752192e-15
[5] 3.5527136788005009e-15 6.8833827526759706e-15
[7] 1.4210854715202004e-14 2.8199664825478976e-14
[9] 5.6843418860808015e-14 1.1346479311669100e-13
[11] 2.2737367544323206e-13 4.5452530628153909e-13
[13] 9.0949470177292824e-13 1.8187673589409314e-12
[15] 3.6379788070917130e-12 7.2757355695785009e-12
[17] 1.4551915228366852e-11 2.9103608412128779e-11
[19] 5.8207660913467407e-11 1.1641509978232989e-10
[21] 6.6493359939245877e-06 2.5720460687386204e-03
[23] 4.8202266324911647e-02 1.7456547460031679e-01
[25] 2.4936031630254019e-01 2.1428293501123058e-01
[27] 1.4106434838399229e-01 8.1018980443629832e-02
[29] 4.3428433624081136e-02 2.2484450838189007e-02
25回が続くのが4回に1回あることになる。
pythonで25回以上と25回ちょうどになるのを計算させてみた。
その結果、
Over 25
6977459029519597/9007199254740992
= 0.7746535667951356
Just 25
2246038053550679/9007199254740992
= 0.24936031612362342
高速化を狙ってCに移植したら
100万回で暴走。
https://egg.5ch.net/test/read.cgi/hosp/1540905566/132
用意されたゼッケンN(=100)枚以下の参加であった。
無作為に抽出したM(=5)人のゼッケン番号の最大値はMmax(=60)であった。
参加人数推定値の期待値とその95%信頼区間を求めよ"
decken <- function(M, Mmax, N, conf.level=0.95){
if(Mmax < M) return(0)
n=Mmax:N
pmf=choose(Mmax-1,M-1)/choose(n,M)
pdf=pmf/sum (pmf)
mean=sum(n*pdf)
upr=n[which(cumsum(pdf)>conf.level)[1]]
lwr=Mmax
c(lower=lwr,mean=mean,upper=upr)
}
decken(M=5,Mmax=60,N=100)
> decken(M=5,Mmax=60,N=100)
lower mean upper
60.0000 71.4885 93.0000
これをシミュレーションで確認したい。
# simulation
M=5 ; Mmax=60 ; N=100
sub <- function(M,Mmax,N){
n=sample(Mmax:N,1) # n : 参加人数n
m.max=max(sample(n,M)) # m.max : n人からM人選んだ最大番号
if(m.max==Mmax) return(n)
}
sim <- function(){
n=sub(M,Mmax,N)
while(is.null(n)){
n=sub(M,Mmax,N) # 最大番号が一致するまで繰り返す
}
return(n)
}
runner=replicate(1e4,sim())
summary(runner) ; hist(runner,freq=F,col="lightblue")
quantile(runner,prob=c(0,0.95))
cat(names(table(runner)[which.max(table(runner))]))
> summary(runner) ; hist(runner,freq=F,col="lightblue")
Min. 1st Qu. Median Mean 3rd Qu. Max.
60.00 63.00 68.00 71.43 77.00 100.00
> quantile(runner,prob=c(0,0.95))
0% 95%
60 93
> cat(names(table(runner)[which.max(table(runner))])) # 最頻値
60
結果は確認できたけど、もっと高速なシミュレーションアルゴリズムはあるだろうか？
覆面算を □/□ ＊ □/□ ＝ □□/□を解くの使ってみた。
a/b * c/d == ef/g (c>a)として
F = function(fun){
n=1:9
ans=NULL
for(a in n){
for(b in n){
for(c in a:9){
for(d in n){
for(e in n){
for(f in n){
for(g in n){
if(fun(a,b,c,d,e,f,g)){
ans=rbind(ans,c(a,b,c,d,e,f,g))}}}}}}}}
return(ans)
}
で虱潰しに判定
F1=function(a,b,c,d,e,f,g){ #全部の組み合わせが等しくないのを確認
(a/b)*(c/d)==(10*e+f)/g &
a!=b & a!=c & a!=d & a!=e & a!=f & a!=g &
b!=c & b!=d & b!=e & b!=f & b!=g & c!=d &
c!=e & c!=f & c!=g & d!=e & d!=f & d!=g & e!=f & e!=g & f!=g
}
F2=function(a,b,c,d,e,f,g){ # anyDuplicatedで重複なしを判定
(a/b)*(c/d)==(10*e+f)/g & !anyDuplicated(c(a,b,c,d,e,f,g))
}
> system.time(F(F1))
user system elapsed
52.56 0.25 53.38
> system.time(F(F2))
user system elapsed
113.78 0.11 115.81
anyDuplicatedでコードは短くなるが速さが犠牲になった。
truehistで出したヒストグラムの情報をテキストファイルとして保存する方法、教えて頂けませんか？
truehist関数はhist関数と違って戻り値を返さないから無理なのでは
MASS::truehist
でソースを覗いたら最後はinvisible()になっていた。
ここを
return(list(breaks=breaks,h=h,nbins=nbins,xlab=xlab))
にしたら
$`breaks`
[1] -0.001 0.999 1.999 2.999 3.999 4.999 5.999 6.999 7.999 8.999 9.999 10.999 11.999
[14] 12.999 13.999
$h
[1] 1
$nbins
[1] 17
$xlab
[1] "x"
で出力されるけど、
どのパラメータがあればヒストグラムが再現できるのかは不勉強にてわからない。
レスありがとうございます。
すると
truehist のソースの最後の
invisible()
を
return(list(breaks=breaks,est=est))
として、breakの中点を横軸、estを縦軸にするとヒストグラムが再現できるわけですね。
に改造したソースをtruehist0として
midpoint <- function(x){ # c(1,2,3,4) -> c(1.5,2.5,3.5)
n=length(x)
mpt=numeric(n-1)
for(i in 1:(n-1)){
mpt[i]=mean(x[i],x[i+1])
}
return(mpt)
}
x=rnorm(10000)
(data=truehist0(x))
with(data, plot(midpoint(breaks),est,type='h',lwd=5,col='cyan'))
で元のヒストグラムが再現できた。
でhistのソースを表示させてみた。
r <- structure(list(breaks = breaks, counts = counts, density = dens,
mids = mids, xname = xname, equidist = equidist), class = "histogram")
if (plot) {
plot(r, freq = freq1, col = col, border = border, angle = angle,
density = density, main = main, xlim = xlim, ylim = ylim,
xlab = xlab, ylab = ylab, axes = axes, labels = labels,
...)
invisible(r)
}
histでのcountsが
truehistではestと呼ばれているようだ。
estimateの略かな？
データの操作が一番むずいんだよ！
ちなみに、フィッティングさせたいのはレブィフライト確率密度関数です。
最小二乗法でもできた気がする
MASSのfitdistとVGAMのlevyを使うとなんとかなるかも。
やったことないけど。
とりあえず、最小二乗法でやってみた。
ガンマ分布の乱数を近似してみた。
dlevy <- function (x,m,c) sqrt(c/2/pi)*exp(-c/2/(x-m))/(x-m)^3/2
set.seed(123)
dat=rgamma(1e3,1) ; hist(dat,freq=F)
x=density(dat)$x ; y=density(dat)$y
lines(x,y)
f<-function(mc){
m=mc[1];c=mc[2]
sum((y-dlevy(x,m,c))^2)
}
(mc=optim(c(0,1),f, method='N')$par)
curve(dlevy(x,mc[1],mc[2]),add=T,col=2)
fitdistは対応できる分布が限定されているから
のソースを改造しないと無理だな。
sink('print.out')
print(MASS::fitdistr)
sink()
でprint.outに出力してみた。
これもoptimを使っているようだが、最小二乗法なのかどうなのかわからなかった。
統計というより待ち行列理論だね
50分が答えになるから合ってそう
レスありがとうございます。
λ=5
μ=6
N=1e6
sum(rpois(N,λ)*rexp(N,μ))/N
> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.833631
ρ=λ/μ
ρ/(1-ρ)
ρ/(1-ρ)*1/μ
> ρ/(1-ρ)*1/μ
[1] 0.8333333
なのでいいのだろうとは思っていたのですが、時系列のシミュレーションは自信がありませんでした。
>550の設定で12分毎に患者が来院したら、待ち時間は全員０だと思うのだが
ρ/(1-ρ) の公式って正しいんだろうかな？
自己レス
公式は定常状態に達したときという前提での計算なんだな。
MMS = function(n, lamda=5,mu=6,s=1){
rho=lamda/mu
sig=0
for(i in 0:s) sig=sig+rho^i/factorial(i)
p0=1/( sig + rho^(s+1)/factorial(s)/(s-rho) )
ifelse(n >= s, rho^n/factorial(s)/s^(n-s)*p0, rho^n/factorial(n)*p0)
}
E=0
for(i in 0:1000) E=E+i*MMS(i)
> E
[1] 5
こういう類の待ち時間の話。
ある医院では、患者が平均１０分間隔でポアソン分布にしたがって訪ねてくることがわかった。
医者は１人であり、１人の患者の診療にかかる時間は平均８分の指数分布であった。
「平均待ち時間」を5分以下にするには同じ診察効率の医師が何人に必要か？
その最小人数で「平均待ち時間」を5分以下に保って診療するには１時間に何人まで受付可能か？
公式に当て嵌めれば解けるのだけど
どうやってシミュレーションすればいいのか思い浮かばない。
コイントスやサイコロだとシミュレーションは容易なんだが。
# ある医院に１時間あたり平均５人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
# 1時間あたりの平均診療人数は６人で、一人あたりの診療時間は指数分布に従うとする。
# 診察までの平均の待ち時間は何時間か？
MM1sim <- function(n=40,lambda=5/60,mu=6/60,seed=FALSE,Print=TRUE){
# service starc clock time(ssct) since 9:00
ssct=numeric(n)
# waiting time(w8)
w8=numeric(n)
# service end clock time(sect)
sect=numeric(n)
# arrival clock time(act)
if(seed) set.seed(1234) ;
act=round(cumsum(rexp(n,lambda)))
# duration of service(ds)
if(seed) set.seed(5678) ;
ds=round(rexp(n,mu))
# simulation assuming service starts at 9:00
head(act) # act : arrival clock time
head(ds) # ds : duration of service
# initial values
ssct[1]=act[1] # 9:15 service start clock time for 1st guest
sect[1]=act[1]+ds[1] # 9:25 sevice end clock time for 1st guest
w8[1]=0
for(i in 2:n){
w8[i]=max(sect[i-1]-act[i],0)
ssct[i]=max(sect[i-1],act[i])
sect[i]=ssct[i]+ds[i]
}
if(Print){
print(summary(w8))
hist(w8,freq=FALSE,col="lightblue",main="")
}
invisible(w8)
}
w8m=replicate(1e3,mean(MM1sim(P=F)))
summary(w8m)
# simulation step by step
#
# act[2] # 9:16 arrival clock time of 2nd
# max(sect[1]-act[2],0) # 9:25-9:16 vs 0 = ?sevice for 1st ends b4 2nd arrival
# w8[2]=max(sect[1]-act[2],0) # 9 min : w8ing time of 2nd
# ssct[2]=max(sect[1],act[2]) # 9:25 vs 9:16 = service start clock time for 2nd
# sect[2]=ssct[2]+ds[2] # 9:25 + 8 = 9:33 service end clock time for 2nd
#
# act[3] # 9:17 arrival clock time of 3rd
# max(sect[2]-act[3],0) # 9:33 - 9:17 vs 0 = ?serivce for 2nd ends b4 3rd arrival?
# w8[3]=max(sect[2]-act[3],0) # 16 min : w8ting time of 3rd
# ssct[3]=max(sect[2],act[3]) # 9:33 vs 9:17 = service start clock time for 3rd
# sect[3]=ssct[3]+ds[3] # 9:33 + 11 = 9:44 service end clock time for 3rd
#
患者来院数４０人程度では定常状態に達しないということみたいだな。
１０万人での待ち時間の分布（理論値５０に近い）
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 5.532 28.619 47.540 67.370 528.511
１００人での待ち時間の分布 （シミュレーションの度にばらつく）
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 6.958 16.564 22.200 31.652 109.204
４０人来院を１万回繰り返した平均の分布
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.304 11.764 19.357 25.449 32.536 177.913
ここに上げておきました。
http://tpcg.io/7psmrQ
結局のところ、待ち時間行列の理論を個人医院に適応しても
定常状態（待ち行列の長さが一定）に達していないと実用性がなさそうだな。
シミュレーションはなんとなく機能していた感じ。
ここに上げておきました。
https://www.tutorialspoint.com/tpcg.php?p=7psmrQ
http://tpcg.io/7psmrQ
個人医院で待ち行列の長さが一定になるほどの受診があるとは考えがたいからこっちの方が現実に即しているんじゃなかろうか
と思うが、解析解はどんなふうになるのか想像もつかん。
https://i.imgur.com/tCOoxF7.png
truehistで軸を対数軸にして表示する方法って分かりますか？
答えたの俺だがソース改造できたのかなぁ。
with(hist(c(rnorm(1e6),rnorm(1e5,5,0.5))),plot(mids,log10(counts),type='h',col=4,lwd=10))
truehistだとソース改造すればいい。
plot(log='y')でもいいな。
例
dat=hist(c(rnorm(1e6),rnorm(1e5,5,0.5)))
attach(dat)
plot(breaks[-1],counts,type='s',log='y',ylim=range(counts))
segments(x0=breaks[-1],y=min(counts),y1=counts)
segments(x0=breaks[1],y=min(counts),y1=counts[1])
heatmap.2でDendrogramつきヒートマップを描きたいのですが、カラムの並びを任意に変えたいです。
Dendrogramの細かい並びは変えないように、大きなクラスタの並びを変えたいです。
例えば、（１，２，３）（４，５，６）（７，８，９）とあるのを
（４，５，６）（７，８，９）（１，２，３）とならび変えるのが目的です。
このとき、4,5,6と7，8，9は近いクラスタを形成します。なので、樹形図は崩れないように書き換えられると思っています。
https://www.biostars.org/p/237067/
上記サイトをみると、as.dendrogramのアウトプットをreorderで並び替えてColvに入れるようですが、
うまくいきませんでした。並びは全く変わっていません。
どなたか教えていただけますか。情報に漏れがありましたらご指摘ください。
環境は以下の通りです。
R 3.4.0
Rstudio 1.0.143
Win10になります。
並べ替えだけなら
x=rbind(c(1,2,3),c(4,5,6),c(7,8,9))
x[c(2,3,1),]
x=list(c(1,2,3),c(4,5,6),c(7,8,9))
x[c(2,3,1)]
Wolframで計算したら確かにそうなった。
https://www.wolframalpha.com/input/?i=integral(x-(a%2Bb)%2F2)%5E2%2F(b-a),from+a+to+b
バスの到着時間が平均１０分の指数分布に従うときにランダムにバス停に行ったときの平均待ち時間は１０分。
バスの到着時間が平均１０分の一様分布に従うときにランダムにバス停に行ったときの平均待ち時間は６分４０秒。
バスがきちんと１０分毎に到着するときはランダムにバス停に行ったときの平均待ち時間は５分。
乱数発生させて公式でのシミュレーション
> d2w8 <- function(x){# w8=E[X2]/2E[X]=(V[X]+E[X]^2)/2E[X]
+ c(mean=mean(x),var=var(x),w8=mean(x^2)/mean(x)/2)
+ }
> N=1e6
> d2w8(rexp(N,1/10)) # exp average:10
mean var w8
10.02477 100.67652 10.03377
> d2w8(runif(N,0,20)) # unif average:10
mean var w8
9.997470 33.325065 6.665408
> d2w8(rep(10,N)) # regular interval 10
mean var w8
10 0 5
プログラムしていたら、これに気づかないのがバグの原因だったｗ
> any(NULL)
[1] FALSE
> all(NULL)
[1] TRUE
> switch (3,
+ x =1,
+ x =2,
+ x =3
+ )
[1] 3
> switch (2,
+ x <- 1,
+ x <- 2,
+ x <- 3
+ )
そのふたつは意味が異なる
それぞれちゃんと命令どおりの挙動
> rm(x)
> x=switch(1,x =1)
> x
[1] 1
> switch(1,x<-1)
> x
[1] 1
その二つは異なるアルゴリズムでたまたま結果が同じになっているだけ
=と<-の違いはRやるなら理解しておいたほうが良い
scopeが違うってことですね。
x^0=1
0^0=1とした方が辻褄が合うことが多いけど
Rのこの仕様には何のメリットがあるんだろ？
> any(NULL)
[1] FALSE
> all(NULL)
[1] TRUE
>
神経質すぎ
平成29年の簡易生命表から
f=c(179,28,19,13,9,7,6,5,5,4,4,4,5,7,9,10,11,12,14,16,18,19,19,20,21,22,23,24,25,27,28,30,
31,34,37,39,41,43,46,51,57,63,70,77,83,91,99,109,119,130,142,155,167,178,190,202,216,
233,249,265,282,302,326,354,387,420,457,502,549,598,648,700,761,836,927,1026,1142,1284,
1455,1651,1862,2089,2341,2625,2934,3264,3598,3923,4233,4508,4740,4893,4973,5007,4999,
4729,4314,3797,3222,2634,2071,1566,1135,788,523,761)
m=c(191,31,21,13,10,8,8,8,7,7,7,8,9,11,14,17,21,26,32,37,42,46,49,50,50,50,49,49,50,52,55,
58,60,63,65,67,71,76,82,89,97,104,112,122,134,148,165,183,203,224,246,268,294,324,357,
391,425,461,502,549,601,659,722,792,872,958,1052,1147,1239,1331,1433,1546,1663,1783,
1905,2026,2167,2333,2532,2750,2973,3195,3414,3630,3827,4000,4133,4200,4194,4104,3916,
3681,3388,3046,2669,2272,1875,1494,1145,841,589,392,246,144,79,71)
LE <-function(ndx,Y,N0=10^5){ # life expectancy
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(m,65)
LE(m,61)
例えばベクトルに値を追加すると
先頭が消えて後尾に新しい値をついかして要素数を一定に保つような
一定要素数以下の平均を求めたいのでそういうのが簡単に実現できる方法あればなおよいです
ほかにもっと簡単な方法アレばお湿気てください
mean(tail(x, 50))
f=function(x,y,n=50) mean(tail(append(x,y),n))
これはおかしい
具体的には
timeStr <- "2018-01-07 01:00"
dateTime <- strptime(timeStr, format="%Y-%m-%d %H:%M") #"POSIXlt" クラスオブジェクトに変換
として時間データに変換したものをx軸としてプロットしたものです
たとえば
plot(x=dateTime, y=1)
abline(v=?????)
として縦線を追加したいのですが
みたいな感じにすれば解決しました
かなり遅れすだけど
formals(cor)
>>617
俺の予想
判定するときにNULLは強制的に論理値に変換される
変換されると logical(0) になる
logical(0) は空のベクトル
からのベクトルが渡されて中身をチェックしていく
アルゴリズムとして早く処理するためには
any()はTRUE探しにいって、一個でもTRUEがみつかればその時点でTRUEをかえす
all()はFALSEを探しに行って、一個でもFALSEがみつかればその時点でFALSEをかえす
もちろん最後までみない
空のベクトルが渡されたのでTRUEもFALSEもみつからない、となると
any()ではFALSEとなり
all()ではTRUEとなる
これで辻褄はあう
ちなみに
any(c())
all(c())
でも同じ結果が出る
1行データの構造の設計からデータフレームの変換までどうすればシンプルに実現できるか助言ください。
たとえば
1 a
2 b
というようなデータを
1-a,2-b
というような一行のテキストデータから初めてテーブル構造にするというような形です
x <- "1-a,2-b"
y <- str_split(x, ",")
と試しにやってみたのですがyが行単位のベクトルになるだけでここからどうデータフレームにすればよいかわかりません
何をしたいのか、まったく理解できない。最低限、N, Mが何なのか説明したらどう？
文字列 1-a,2-b を 数値とテキストのデータフレームにしたいという意味と解した。
x="1-a,2-b"
y=strsplit(x,",")
z=unlist(y)
w=NULL
for(i in 1:length(z)){
w=rbind(w,unlist(strsplit(z[i],"-")))
}
data.frame(NUM=as.numeric(w[,1]),TEXT=w[,2])
> x="1-a,2-b"
> y=strsplit(x,",")
> z=unlist(y)
> w=NULL
> for(i in 1:length(z)){
+ w=rbind(w,unlist(strsplit(z[i],"-")))
+ }
> data.frame(NUM=as.numeric(w[,1]),TEXT=w[,2])
NUM TEXT
1 1 a
2 2 b
>
NMは任意の数字です
したいことは
データを一行に記録してそれを
テーブル構造にすることです。
これだけです。
一旦ファイルに保存してread.csvなどにすればよいのでしょうが
いちおう直接テキストデータをコピペしてからということにしたいのです。
なぜこんなことをするかと言うと
TamperMonkeyという拡張機能でJavaScriptで使ってウェブ上のデータを収集しているのですが、
localStrageというクッキーの拡張版のような機能をつかってデータをクライアントに保存するときに基本的にテキストデータベースでの保存になので
一行のデータとして後ろにどんどんデータを追加していくのが一番単純な処理に成るからです。
そのlocalStorageに保存されたデータはコピペで取り出すしかないので
それをRで処理する時に一行のテキストデータから始めないといけないのです。
, を行のデリミタ、- を列のデリミタとして、一行の文字列をデータフレームにするというのであれば、
s に文字列が入っているとして、
r <- unlist(strsplit(s, ","))
d <- lapply(r, function(x) unlist(strsplit(x, "-")))
as.data.frame(do.call(rbind, d))
列数が行によって異なるときにどうなるかは知らん。
内部動作の考証ありがとうございます。
未だにこれは理解できません
> logical(NULL)
Error in logical(NULL) : invalid 'length' argument
> logical(0)
logical(0)
> logical(1)
[1] FALSE
> txt <- "1-a,2-b,3-c"
> read.table(text = gsub(',', '¥n', txt), sep = '-')
V1 V2
1 1 a
2 2 b
3 3 c
こんな感じか？
自分の解釈では、
NULL は「無」なのでエラー（引数はベクトルの要素数を必ず与えなければいけない）
0のとき、0個という指定なので、0個の要素をもつベクトル
1のとき、1個という指定なので、1個の要素を持つベクトル（規定値はFALSE）
ちなみに、
> logical(2)
[1] FALSE FALSE
2のとき、2個という指定なので、２個の要素を持つベクトル（規定値はFALSE）
解説ありがとうございました。
>>637
ファイル入れなくてもできるんですね
一巡回って大化2でいいよ
あとは干支みたいに回せばいい
レス付くと思わなかった。。。
⇣の結果をすべて即答できるようになれば合格だと思います…
any(c(TRUE,NA))
any(c(FALSE,NA))
any(c(FALSE,NA,TRUE))
any(c(NA, TRUE))
any(c(NA,FALSE,TRUE))
all(c(TRUE,NA))
all(c(FALSE,NA))
all(c(FALSE,NA,TRUE))
all(c(NA, TRUE))
all(c(NA,FALSE,TRUE))
TRUE | NA
FALSE | NA
FALSE | NA | TRUE
NA | TRUE
NA | FALSE
NA | FALSE | TRUE
TRUE & NA
FALSE & NA
FALSE & NA & TRUE
NA & TRUE
NA & FALSE
NA & FALSE & TRUE
例えば、下記のコードでは線と点で別々の凡例になります。
線と点のスタイルを合わせて１つの凡例にしたいのですがどうすればいいでしょうか。
geom_line(mapping=aes(colour=Conditions),alpha=0.6)
+ geom_point(mapping=aes(shape=Conditions,colour=Conditions),alpha=0.8)
+ scale_shape_manual(values = 1:3)
このコードではなくもっと後のthemeで凡例のスタイルを変えるときに余計なことをしていました。
mean(a)
#> [1] "a"
これ本に書いてたんだけど
君たちこれ実行して"a"がでてくる?
自分とこでじっこうしても
> mean(a)
[1] NA
警告メッセージ:
mean.default(a) で: 引数は数値でも論理値でもありません。NA 値を返します
となるんやけど?
いみわからん。
http://adv-r.had.co.nz/OO-essentials.html
でもJSでOOかじった程度なので全然意味分からん
なにやってんのこれ?
オブジェクト指向をちゃんと勉強してくれ。
横からだが、勉強になった。
クラスの自作とか考えたことがなかったけど、今度、俺様クラスを作ってみよう
知ってるかもだけどいちいちmappingは書かなくてもいい
scale_shapeも値が1:3ならいらない
ggplot(data, aes(shape=Conditions, col=Conditions))+
geom_line(alpha=0.6)+
geom_point(alpha=0.8)
https://logtube.jp/variety/28439
【衝撃】ヒカキンの年収・月収を暴露！広告収入が15億円超え！？
https://nicotubers.com/yutuber/hikakin-nensyu-gessyu/
HIKAKIN(ヒカキン)の年収が14億円！？トップYouTuberになるまでの道のりは？
https://youtuberhyouron.com/hikakinnensyu/
ヒカキンの月収は１億円！読唇術でダウンタウンなうの坂上忍を検証！
https://mitarashi-highland.com/blog/fun/hikakin
なぜか観てしまう！！サバイバル系youtuberまとめ
http://tokyohitori.hatenablog.com/entry/2016/10/01/102830
あのPewDiePieがついに、初心YouTuber向けに「視聴回数」「チャンネル登録者数」を増やすコツを公開！
http://naototube.com/2017/08/14/for-new-youtubers/
27歳で年収8億円 女性ユーチューバー「リリー・シン」の生き方
https://headlines.yahoo.co.jp/article?a=20170802-00017174-forbes-bus_all
1年で何十億円も稼ぐ高収入ユーチューバー世界ランキングトップ10
https://gigazine.net/news/20151016-highest-paid-youtuber-2015/
おもちゃのレビューで年間12億円！ 今、話題のYouTuberは6歳の男の子
https://www.businessinsider.jp/post-108355
彼女はいかにして750万人のファンがいるYouTubeスターとなったのか？
https://www.businessinsider.jp/post-242
1億円稼ぐ9歳のYouTuberがすごすぎる……アメリカで話題のEvanTubeHD
https://weekly.ascii.jp/elem/000/000/305/305548/
世界で最も稼ぐユーチューバー、2連覇の首位は年収17億円
https://forbesjapan.com/articles/detail/14474
ヒカルの収入が日収80万、月収2400万、年収3億と判明www
https://matomenewsxx.com/hikaru-income-8181.html
はじめしゃちょーの年収は6億？2017年は30億突破か？
https://2xmlabs.com/archives/1873
2月1日とか聞いたような
可読性を保つために省略しないほうがいい
プログラミングの基本
特にmappingなんかは書くまでもない。
マウント取りもいきすぎると滑稽だわ
お前のほうが滑稽
マウントってるのお前
イキってるだけで中身もないし
滑稽すぎる
df［-1:-grep("start",df[,1])］
今こんな感じだけどパイプで繋げない。
slice使えばいいのか、grep以外にいい関数があるのか？
startが必ず一つしかないとか、restartとかがあったときにそれを検出しても良いならそれでもよいが、普通は
df[-(1:which(df[, 1] == "start")[1]), ]
すると思う。
経験的にはgrepはできるだけ避けた方がいい。
スレ違いが明らかなんですが、ここで聞かせて下さい。
numpyスレ立てていいですか?
とんでもなくスレ違い
さっさと立てるね
なるほど、ありがとうございます。
そしてカッコでくるんでからマイナス付ける書き方もあるのか。勉強になります。
可読性について勉強できてよかったな
現状col<-apply(2,function(x)｛all(is.na(x))｝
で要らない列定義してからdf[,!col]としてるけどほんとはパイプの中に入れて処理したい。
自己解決しました。select_ifでいけそうです。
板違いだと思うので、やめておいた方が良いのでは
ノーマルモードのときに日本語入力してしまうとバグみたいになるんだが
あれどうにかならんの?
なんとかしてくれんかなぁ
IMEが無効になるバグが直ってる感じがする。
まだ十分使い込んでないからだけかも。
シロンボヒトモドキゴキブリニホンザルの自由は偽自由と詐欺広告人を殺す自由ヒトモドキゴキブリシロンボアメ公はニホンザルゴキブリと自殺せよ？
あほな反抗期の中学生がつかってるイメージしかわかんのだが。
いずれにしろ広辞苑に載ってるから調べればいい
http://mevius.5ch.net/test/read.cgi/seiji/1538271437/l50
まもなく世界経済破綻、みなさんの財産は紙くずになります。
方言に対して広辞苑を持ち出すのは的外れ。
大阪弁辞典によると「意気る」または「粋る」なので、漢字も外れ
どっちも「熱る」の当て字やん
「熱る」って初めて聞いた
これは普通に使われているの？
方言？
生まれてこのかた、東京なんだけど聞いたこと見たことないんだが
単なる私の不勉強かな？
いき・る【熱る】
１ あつくなる。ほてる。むしむしする。
２ 激しく怒る。
https://kotobank.jp/word/%E7%86%B1%E3%83%BB%E7%86%85-2005439
いき・る【熱る】
① あつくなる。ほてる。むしむしする。
② 息づかいを荒くして怒る。相手と争おうとしていきまく。言いたてる。
③ 調子に乗って勢いこむ。元気づく。
ヒトモドキ奇形なりすましシロンボアメ公はイスラム国に首切られろヒトモドキフィフィラクダイスラム国の売女とともに自殺しろテロゴキブリ民族アメ公
キモ豚山田太郎は自民入党のマイノリティ下僕化しキモオタヒトモドキ障害者の存在価値は0に達したんだよ表現戦士の性犯罪者害虫
このキモオタヒトモドキネトウヨをガス室で抹殺してえよな？w
奇形ネトウヨヒトモドキゴミ藪部落の顔気持ち悪いから自殺しろ糞犬hk
ヒトモドキ産経便所ゴキブリは下着ドロの常習犯変態暴論ゴミ便所自殺しなさい今すぐ
インフルエンザの新しい治療薬「ゾフルーザ」を投与されたＡ香港型のインフルエンザ患者30人を調べたところ、70％余りに当たる22人から、この薬が効きにくい耐性ウイルスが検出されたことが国立感染症研究所の調査で分かりました。
調査件数は多くないものの、専門家は現在のような使用を続けると、耐性ウイルスが広がるおそれがあるとして使用基準を見直すべきだと指摘しています。
耐性化率が50%以上である確率は
pbeta(0.5,1+22,1+8,lower=F)
[1] 0.9946631でいいかな？
[1] 0.9838752
> binom.test(22,30,conf=p)
Exact binomial test
data: 22 and 30
number of successes = 22, number of trials = 30, p-value = 0.01612
alternative hypothesis: true probability of success is not equal to 0.5
98.38752 percent confidence interval:
0.5000000 0.8994036
sample estimates:
probability of success
0.7333333
google料金高過ぎて鞍替えなんですが
使ってるQtが古いんや。1.2でだいぶましになった感じだけど入力途中の文字が残る妙な挙動が時々でる。
書籍は出版される頃には情報が古くなってるから。自分はネットの情報とヘルプだけで十分だな。
ダウンロードしようとしても32ビット版見当たらないし、Windows7,10用の64ビット版セットアッププログラムを起動しようすると32ビットへの選択とかはなくて、単にエラーになっちゃうし
誰か教えてください
低脳中卒万引きヒトモドキネトウヨ猿ヒトモドキをさしころせ
古いのはここにあるけど今後のことをかんがえたらOSを64bit版に変えた方がいいよ。
https://support.rstudio.com/hc/en-us/articles/206569407-Older-Versions-of-RStudio
学生ですか？
仕事あれば教えてほしい
Rを仕事にしてはいないけど、仕事の一環でR使ってる
信頼性、入手しやすさ、扱いやすさのバランス的に自分の環境ではR一択なんだよな
ちなみにPythonはどう？
AI目当てでこれから使おうと思って勉強中
今の職場で求められるのは主に統計用途や作図なんでRは手放せないかな
ゴキブリネトウヨヒトモドキ奇形売春婦肉便器なつこババア滅多刺しにして解体しろ
なるほど。
自分も機械学習やってるけど、R の方がやりやすいよ。ディープラーニングならPythonかもしれないけど
機械学習よりは統計寄りだとRは仕事でもスタンダード、特にマニアックな統計解析だとRにしかパッケージがない
まあ使い慣れてるならRで機械学習やっても全然構わないけど
サンクス
まあメモリが2Gなので
少し前はRなんてだんだん廃れるみたいな雰囲気あったけど今はむしろ盛り返してデファクトな何時になってるからこの先10年は行きそうだなあ
機械学習特需でpythonの躍進はすごかったけど
データファイル名を変えてフォルダには変更名で保存されてるのに、読み込むともとの名前のままなんです。どうすれば変更後の名前になるのでしょう？
サンクス
当面は32ビット版で頑張る
全く状況がわからん。
主語や目的語、必要な修飾語を省略せずに、分かるように説明しないと誰も助言できない。
なんなんだろう
それはsave, loadよ話か？
ふうL@Fu_L12345654321
学コン1傑いただきました！
とても嬉しいです！
ｈｔｔｐｓ://pbs.twimg.com/media/D-IuUuqVUAALnAB.jpg
ｈｔｔｐｓ://twitter.com/Fu_L12345654321/status/1144528199654633477
https://twitter.com/5chan_nel (5ch newer account)
> 【備えあれば】防災用品 非常食スレ121【憂いなし】
> https://mao.5ch.net/test/read.cgi/eq/1562388281/
今から始めるならどっちがオススメ？
文字化け？
R本体だとグラフィックディバイスに違いはあるけど、実質的に差はない。
でもRは本体だけで使うことは稀で、多種多様なGUIがあるから、
そこでOSの違いによる差が生じる。
# サイコロ
# 正６面体のサイコロがある．４面は青色、２面は赤色である．
# このサイコロを合計20回振るとき、最も起こりそうな順番はどれか？
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
sim <- function(){
a=sample(0:1,20, replace=TRUE, prob=c(4,2))
b=as.character(a)
c=paste(b,collapse="")
s1=paste(c(1,0,1,1,1),collapse="")
s2=paste(c(0,1,0,1,1,1),collapse="")
s3=paste(c(0,1,1,1,1,1),collapse="")
res=c(grepl(s1,c),grepl(s2,c), grepl(s3,c))
return(res)
}
k=1e6
re=replicate(k,sim())
mean(re[1,])
mean(re[2,])
mean(re[3,])
結果は、直感とおり、１が再頻
> mean(re[1,])
[1] 0.124672
> mean(re[2,])
[1] 0.080873
> mean(re[3,])
[1] 0.040564
>
grep使わない方法ってあるかな？
ワイの直感的解法
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
だが、以下でも確率同じ。何となく
# １．赤 赤 赤 赤 青 ★
# ２．赤 赤 赤 赤 青 青
# ３．赤 赤 赤 赤 青 赤
★は赤でも青でもどっちでもOK
P(# １) >P(# ２) >P(# ３) ∵青出やすい
R言語等でシミュレーションされ、
自身の確率直感が正しいのを
確認できるとは、素晴らしい。
# １．赤 青 赤 赤 赤
# １'．赤 赤 赤 赤 青 とすると、
# １と# １'は、直感で同じ確率と
思ってたが間違えのようだ。
当方のシミュレーションで、
# １は、0.1248
# １'は、0.1271
となった。微妙だけど、多分だ。
やっぱり確率計算をコンピュータで
モンテカルロシミュレーションのは
素晴らしい。
【grepは未使用の糞真面目な方法】
# １．赤 青 赤 赤 赤
についての確率、ほぼ厳密解を得た
# １は、0.124774… だと思う
計算は、モンテカルロ法でない方法
でプログラム、計算した。
で、grepは使用してない。
ちなみに計算誤差は、ほぼ皆無なハズ
ソースコードイメージ
p01 = (1/3)^4*(2/3)
p02 = p01
p03 = p01 * (1 - p01)
p04 = p01 * (1 - p01 - p02)
p05 = p01 * (1 - p01 - p02 - p03)
…
p16 = p01 * (1 - p01 - p02 - p03 … - p14)
とし、
p01~p16 の合計を算出したところ、
0.124774… となった
レスありがとうございます。
私の直感
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
# １．を６個に書き換えて #２．と並べると
# １．★ 赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
★は赤でも青でもどっちでもOKだから#１．の方が起こりやすい
# ２と# ３を比べると
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
３個めでは青の方がでやすいので
ら#２．の方が起こりやすい
よって、P(# １) >P(# ２) >P(# ３)
>>
1回3.6%で激レアが出るガチャを10回回した確率って
36%なのでしょうか？
それとも0.964*0.964*0.964(略 0.964を10回電卓にかけた数なのでしょうか？
教えてください。
<<
百万回シミュレーション
p=3.6/100
N=10
sim <- function() any(rbinom(N,1,p)==1)
mean(replicate(1e6,sim()))
> mean(replicate(1e6,sim()))
[1] 0.306628
gacha=c(rep(1,36),rep(0,1000-36))
sim2 <-function() any(sample(gacha,10,replace=TRUE)==1)
mean(replicate(1e6,sim2()))
> mean(replicate(1e6,sim2()))
[1] 0.306904
何故か0.36にならない、どうしてだろ？
gacha=c(rep(1,36),rep(0,1000-36))
sim3 <- function() any(replicate(10,sample(gacha,1))==1)
mean(replicate(1e5,sim3()))
> mean(replicate(1e5,sim3()))
[1] 0.30691
これも0.307弱だな。
何が悪いんだろ？ 俺の頭かな？？
(1 - 0.036) ^ 10 ≒ 0.6930592
したがって、
3.6%のガチャを10回まわして1回以上当たる確率は
1 - 0.6930592 = 0.3069408
ありがとう、
やっぱ、
これだったｗ
何が悪いんだろ？ 俺の頭かな？？
シミュレーションにバグはなかったｗ
事実、結果こそすべて
激レアは、レアだから1個とみなす。
故に、母集団の個数は
1÷0.036 = 27.777… きっと28個だ
手順1) 28枚のカードがある
手順2) 重複しない1~28の番号を振る
手順3) 28枚のカードをシャッフル
手順4) 1~10枚目のどれか1となる確率
と絶対同じハズ、だから、
P(1枚目で当) =27P27 ÷ 28P28 = 1/28
P(2枚目で当) も同様に1/28
P(3枚目で当) も同様に1/28
…
P(10枚目で当) も同様に1/28
で此等10個の事象は背反事象だから、
P = 10/28 = 0.357
∵有効数字3桁と勝手にしちゃう
ところてガチャってゲーム何か
よくわかんないけど、計算しちゃった
そのガチャ3.6%の件、
超幾何分布なのか。二項分布なのか。
確率の小さいとか、母集団が小さい とかだと無視できないと思われる。
確率統計はギャンブル派生数学ぢゃ。
生半可な知識ではカモにされる。
現代の若者たちは、数学特に
確率統計はじめとするギャンブル
の能力が特段に欠けており、
R言語等のプログラミング教育で
ギャンブルゲームを学習すべきだ。
奇麗事の学問だけの今日の数学ぢゃ
カモにされるだけ。
健全な娯楽として賭博系確率統計学
をC R Java Pyson Javascript BASIC
の何れかを学校で学習すべきだ。
かつて、昭和の駄菓子屋によくある
ガチャガチャすなわちカプセルトイ
そのような健全的なギャンブルが
超幾何分布とか二項分布の理解に
役立つのだ。
限られたお小遣の10円玉数枚で、
如何にレアアイテムをGetするかを
子供らは、思考するからだ。
「残り物には福がある」は
確率統計的には正しいのか
子供同士で文学的に議論したものだ。
さて、今のガチャは恐らくは、
デジタルの媒体のスマホゲームだ。
二項分布でよいだろう。
R言語には、
二項分布の密度関数、それの累積関数
はモチロン、それに従う乱数生成を
提供してるようだ。
ゲームの仕組みが複雑化する今、
R言語等の乱数生成プログラミングで
これからのデジタル化ギャンブル社会
で、お金より大切な激レアをドンドン
無限にゲットできる人材の増加を
期待できる可能性を秘めている
同じカプセルトイが10台ある。カプセル取り出し後は補充されない。
アタリを１個でも手に入れる確率は
１台から10個取り出す場合(G10)と
１台から１個を１０台で取り出す場合(G01)
ではどちらが高いか？
そのシミュレーション
rm(list=ls())
N=1000 ; K=36 ;n=10 # アタリ3.6％
g=rep(c(1,0),c(K,N-K))
G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元（超幾何分布）
G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
mean(replicate(1e6,G10()))
mean(replicate(1e6,G01()))
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.307745
[1] 0.3081121
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.307295
[1] 0.3069408
１台から１０個取り出した方がいいみたい。
どの台も全部で、1000カプセルで、
どの台も当りが、36カプセルだと
G10(1台で10個)取り出す方が
僅かですが、確率良さそうですね。
シミュレーションも理論値も
同様な結論のようであり、
G10(1台で10個)取り出す方が僅かに
有利が分かり何か楽しかったです。
仮に、もしどの台も
250カプセル中当たり9カプセルなら、さらにG10戦略がG01戦略より有利な
感触を掴めました。
1000個だとシミュレーションの差が微妙で再現性が不安だったけど
250個にすると、差が明らかにつきますね。
> N=250 ; K=9 ;n=10 # アタリ3.6％
> g=rep(c(1,0),c(K,N-K))
>
> G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元（超幾何分布）
> G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
>
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.31082
[1] 0.3117069
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.306782
[1] 0.3069408
これと rbinomを使うと sample関数を使わなくても同じことができた。
N=250 ; K=9 ;n=10 # アタリ3.6％
mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n) # 非復元（超幾何分布）
mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n # 復元(二項分布)
結果
> mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.311581
[1] 0.3117069
> mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n
[1] 0.307213
[1] 0.3069408
1台から順次１０個取り出すのをイメージすれば残り物には福があるとも言えなくもないな。
>
論の始まりだと言われている。当時の有名な賭博師メレがパスカルに以下のような問題を持ち
込み、その問題についてがその手紙のやりとりの中で論じられているそうである。
甲乙二人がおのおの32ピストル（当時のお金の単位）の金を賭けて勝負したとする。
そしてどちらかが先に３点を得たものを勝ちとし、勝った方がかけ金の総額64ピストルをもら
えるとする。ところが甲が２点、乙が１点を得たとき、勝負が中止になってしまった。
このとき、二人のかけ金の総額64ピストルを甲と乙にどのように分配すればよいだろうか。
ただし二人の力は互角で、勝つ確率はそれぞれ1/2ずつだとする。"
# 先にw(=3)点を得たものを勝ち、甲がA(=2)点、乙がB(=1)点を得たとき、勝負が中止
# 甲の確率を求める
gambling <- function(A=2,B=1,w=3,k=1e5){ # k : how many simulate
sim <- function(){
while(A < w & B < w){
g = rbinom(1,1,p=0.5)
if(g==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}
> gambling(2,1,3)*64 #64ピストル配分
[1] 48.13056
> gambling(2,0,4) # 日本シリーズ
[1] 0.80945
で、
自己流の直感でラフな計算をしてみた。
甲1点、乙1点 ⇒1:1按分 ∵自明 ★
甲3点、乙1点 ⇒1:0按分 ∵自明 ☆
安直に★と☆の平均を考えると
甲2点、乙1点 ⇒1:0.5按分 ∵安直
で、総額で64ピストルだから、
甲に、64*1/1.5 = 42.67 ピストル
乙に、64*0.5/1.5 = 21.33 ピストル
を配分すれば、思ってしまった。
ラフな計算だと、少々不公平なことが
起こるようだ。
R言語等のシミュレーションは、
超天才のワィの直感より、
さらに厳密な確率計算ができそうだ。
最後はAが得点して勝負が決定するのを忘れずにAが勝者になるまでのBの点数で場合分け
NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B：現在の点数 ,p:甲の勝率
ans=0
for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数
ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A)
}
return(ans)
}
> NS(3,2,1,0.5)*64 #64ピストルの分配
[1] 48
> NS(4,1,0,0.5) #日本シリーズで先勝したほうが優勝する確率
[1] 0.65625
NSは日本シリーズの頭文字から命名、No Skinではないｗｗ
ワィの脳内は活性化される だから、
ドンドン、妄想全開となるぅぅぅ。
という訳で、妄想内容を記載する
17世紀の地球に行き、胴元になって
賭博ビジネスでガッポリ儲けるゼェ。
17世紀の地球で、胴元になって、
コイントス賭博事業をするゼェ。
で、システムは、2回コイントスし、
表の回数が0回⇒参加者が32万円GET
表の回数が1回⇒胴元に、32万円PAY
表の回数が2回⇒参加者が32万円GET
参加料は、1万円
まんまと、参加者が沢山 1000人いた
ある参加者は、以下の様に怪説をした
組み合わせは、3通りだ。
1通り目 表の回数が0回
2通り目 表の回数が1回
3通り目 表の回数が2回
だから、参加者の勝つ確率は2/3だ
で、モチロン、胴元ワィの商売は成功
およそ、1000万円儲かっちゃた。
しろ@huwa_cororon 11月27日
苦節6ヶ月、初満点&一等賞です！
https://twitter.com/huwa_cororon/status/1199593474128896000
https://twitter.com/5chan_nel (5ch newer account)
日本シリーズは先に4勝したチームが優勝。
勝率はそれまでの（引き分けを除いた）通算勝率に従うとする。
シリーズ開始前の通算成績はA:2勝、B:4勝であった。
今シリーズでAが先勝（第一試合に勝利）した。
この時点でA,Bのどちらに賭ける方が有利か？
https://rio2016.5ch.net/test/read.cgi/math/1567866548/845
>>
845 名前：１３２人目の素数さん[] 投稿日：2019/12/21(土) 05:52:50.74 ID:wum1jR1j
...
10進数で表された5/9を3進数になおせという問題があって、答は0.12だったのですが、
...
<<
こんなのを作ってみた。
# 小数点付きの数numをN進法で表示する
rm(list=ls())
dec2basen <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納
x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
if(N<=10){ # Nが10以下は数値として表示
cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='')
}
else{ # Nが1１以上は整数部分と小数部分を数列で表示
print(list(int=r,decimal=a))
}
invisible()
}
> dec2basen(5/9,3)
0.120000000000000
> dec2basen(3.14159265359,7)
3.06636514320
> dec2basen(3.14159265359,8)
3.11037552421
> dec2basen(3.14159265359,9)
3.12418812407
> dec2basen(3.14159265359,14)
$int
[1] 3
$decimal
[1] 1 13 10 7 5 12 13 10 8 1 4
バグの原因を追求したら、Rの仕様のせいみたいだな。
> (1.2-1)*5==1
[1] FALSE
> (1.2-1)*5>1
[1] FALSE
> (1.2-1)*5<1
[1] TRUE
これがfloor関数が誤値を返してくる原因だった。
> floor(1)
[1] 1
> floor((1.2-1)*5)
[1] 0
数値演算についてもっと勉強しな。
たいていの言語でそうなる。
IPython 6.5.0 -- An enhanced Interactive Python.
C:\bin\Anaconda3\lib\site-packages\ipykernel\parentpoller.py:116: UserWarning: Parent poll failed. If the frontend dies,
the kernel may be left running. Please let us know
about your system (bitness, Python, etc.) at
ipython-dev@scipy.org
ipython-dev@scipy.org""")
(1.2-1)*5==1
Out[1]: False
(1.2-1)*5>1
Out[2]: False
(1.2-1)*5<1
Out[3]: True
(1.2-1)*5
Out[4]: 0.9999999999999998
> 0.728*5-3.64
[1] -4.440892e-16
> round(0.728*5-3.64,5)
[1] 0
んで、デバッグ版
# 小数点付きの数numをN進法で表示する(62進法まで対応 0-9,a-z,A-Z)
rm(list=ls())
dec62 <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示数列を格納
k=max(nchar(num)-nchar(floor(num))-1,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
x=round(num-floor(num),k) # e.g. 7.28-floor(7.28)-0.28 != 0に対応
for(i in 1:k){
y=round(x*N,k) # e.g. 0.728*5-3.64 !=0 に対応
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters,LETTERS)[1:N]
if(N<=62){ # Nが62以下は数値として表示
cat(paste(fig[b$integer+1],sep='',collapse=''),
'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが63以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b) # b$num:検証用
}
Haskellでもやってみた
Prelude> (1.2-1)*5
0.9999999999999998
Prelude> 0.72*5-3.6
-4.440892098500626e-16
レスありがとうございました。
アドバイスに従ってちょっと勉強してみました。
他の言語に移植しても無駄とわかって助かりました。
こうなる理由が理解できました。
(1+1/10-1)*10==1
[1] FALSE
> (1-1+1/10)*10==1
[1] TRUE
これをwhile使って高速化するにはどうすればいいだろう？
俺の頭では思いつかない。
perm <- function (n) {
v=1:n
sub <- function(n, v) {
if (n == 1)
matrix(v, 1, 1)
else {
x = NULL
for (i in 1:n) x =rbind(x, cbind(v[i], sub(n -1,v[-i])))
x
}
}
sub(n, v)
}
> perm(4)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 1 2 4 3
[3,] 1 3 2 4
[4,] 1 3 4 2
[5,] 1 4 2 3
[6,] 1 4 3 2
[7,] 2 1 3 4
[8,] 2 1 4 3
[9,] 2 3 1 4
[10,] 2 3 4 1
[11,] 2 4 1 3
[12,] 2 4 3 1
[13,] 3 1 2 4
[14,] 3 1 4 2
[15,] 3 2 1 4
[16,] 3 2 4 1
[17,] 3 4 1 2
[18,] 3 4 2 1
[19,] 4 1 2 3
[20,] 4 1 3 2
[21,] 4 2 1 3
[22,] 4 2 3 1
[23,] 4 3 1 2
[24,] 4 3 2 1
>
自己解決 e1071というパッケージに再帰呼び出しなしのpermutationがありました。
library(e1071)
permu <- function (n) {
z <- matrix(1)
for (i in 2:n) {
x <- cbind(z, i)
a <- c(1:i, 1:(i - 1))
z <- matrix(0, ncol = ncol(x), nrow = i * nrow(x))
z[1:nrow(x), ] <- x
for (j in 2:i - 1) {
z[j * nrow(x) + 1:nrow(x), ] <- x[, a[1:i + j]]
}
}
return(z)
}
Cで
#include <stdio.h>
int main(int argc, char* argv[])
{
while(argc--)
printf("%s\n", *argv++);
return 0;
}
をRでやるには ... を使って、c(...)もしくはlist(...)で引き出せばいいんだな。
検索してもなかなかわからなかった。
main <- function(...){
argv=c(...)
for(i in 1:length(argv)) cat(argv[i],'\n')
}
library(Rmpfr)
mpfr(factorial(50),5000)
one = mpfr(1, 5000)
factorial(one*50)
> mpfr(factorial(50),5000)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713018969967457666435945132957882063457991132016803840
> one = mpfr(1, 5000)
> factorial(one*50)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713378043612608166064768844377641568960512000000000000
何故か、oneを高精度に設定しておくと、正しい計算をしてくれた。
> round(2.5)
[1] 2
> round(1.5)
[1] 2
は丸め誤差が減るという説明だったので実験してみた。
# 四捨五入 vs ＩＥＥＥ with FUN(=mean) for x(=c(-0.9,-0.8,...,0.8,0.9))
comp <- function(n=10,FUN=mean,x=sample((-9:9)/10,n,rep=T),print=T){
X=FUN(x)
dif = abs(X-FUN(round(x))) - abs(X-FUN(f45(x))) # round後の実行と四捨五入後の実行の差
if(dif!=0 & print) cat(dif<0,' : ',sort(x),'\n') # 差があれば表示 round後が小さければTRUE
return(dif) # 差を返す dif<0:round優位 dif>0:四捨五入優位
}
comp()
k=1e5
# mean
re=replicate(k,comp(print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
# squared sum
comp(FUN=function(x) sum(x^2))
re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
> # 平均 mean
> re=replicate(k,comp(print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ 四捨五入 引き分け
0.25235 0.16199 0.58566
> # 平方和 squared sum
> re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ 四捨五入 引き分け
0.40519 0.01413 0.58068
平均だと小差だが、統計計算の内部処理で使う平方和だと大差がついたので、体感できた。
プラスで14日隔離しないといけない。それが今の船内の状況という。
両隣のどちらかが感染したら14日延長、どの部屋も1日で感染する確率pは1%
部屋の配置は長方形でどの部屋にも両隣があるとする。
発症するか、隔離期間が終われば下船できる。
全員定員１の個室として客と乗務員を合わせた人数nは3000人。
クルーズ船から全員下船できる日数の期待値は？"
p=0.001 # 1日の感染確率
n=3000 # 収容者人数
d14=14 # 隔離日数
Q0=rep(d14,n) # 必要隔離日数の配列の初期値
#
ODA <- function(X){ # One Day After 必要隔離日数の配列 x から１日後の配列を返す
Q=X # Xを作業用Qに代入
S=which(Q!=0) # 未感染者susceptibleのindex
s=length(S) # 未感染者人数
m=rbinom(1,s,p) # その日の感染数
if(m==0){ # 感染者0なら未感染者の
Q[S]=Q[S]-1 # 必要隔離日数を１日減らして返す
invisible(Q)
}else{
I.idx=sample(S,m) # 未感染者Sからm人が感染、感染者のindex (Infected.idx)
id2ext <- function(id){ # 感染者idからの両隣のidを返す
plus =ifelse(id==n,1,id+1) # 最大番号なら1を返す,それ以外は1増やす
minus=ifelse(id==1,n,id-1) # １なら最大番号を返す,それ以外は１減らす
unique(c(plus,minus)) # 重複を除く
}
E.idx=as.vector(sapply(I.idx,id2ext)) # 感染者の両隣のidを列挙し
E.idx=unique(E.idx) # 重複を除く Extended quarantine.idx
Q[E.idx]=d14 # 両隣は必要隔離日数をd14にリセット
Q[I.idx]=0 # 感染者は必要隔離日数=0 (リセット者に重複していたら０で上書き)
IE.idx=unique(c(I.idx,E.idx)) # 感染もしくは感染者の隣のindex
minus1ifnotzero <- function(x) ifelse(x==0,0,x-1) # ゼロでなければ１を減じて返す
Q[-IE.idx]=minus1ifnotzero(Q[-IE.idx]) # 感染もなく感染者の隣りでもない人の隔離日数を１日減らす
invisible(Q)
}
}
sim <- function(){
AZ=FALSE # All Zero 隔離日数がすべて０か？
i=0 # カウンター
Q=Q0 # 隔離日数配列初期値
while(!AZ){ # 隔離日数がすべて0になるまで
Q=ODA(Q) # One Day After関数を実行
AZ=all(Q==0) # All zero？
i=i+1 # カウンター
}
return(i) # 何日かかったかを返す
}
k=1e4 # k回のシミュレーションの
RE=replicate(k,sim())
mean(RE) # 平均値＝期待値
低すぎぃ
一日の感染確率。
3000人いて1日3人くらい報告されたから0.1%/日くらいかな
コードにまだバグがあるかも。
csvそのものをutf-8にしたら今度は不正なマルチバイトがあるってエラー
どうすればいいか教えて
どんなデータ？
ヘッダーだけじゃなくてセル内も日本語？
エスパーじゃないから、まず確認だけど
・Rを使ってる環境(OS)とRのバージョンは？
・CSVのコード変換に使ってるソフトは？
・CSVを読み込んでる関数と指定している引数は？
後で遊んでみよう
エラーがでた。R 3.6.3にアップデートしたけど
> library('remotes')
> remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
Downloading GitHub repo GuangchuangYu/nCov2019@master
Skipping 1 packages not available: BiocStyle
Installing 5 packages: downloader, cowplot, maps, magick, BiocStyle
Error: Failed to install 'nCov2019' from GitHub:
(converted from warning) package ‘BiocStyle’ is not available (for R version 3.6.3)
エラーと通り、BiocStyleパッケージがないので別途、bioconductorからインストールすればいい。
https://bioconductor.org/packages/release/bioc/html/BiocStyle.html
ありがとうございます。
という記述があるんだけど 関数 nls による非線形回帰分析 で作成しているのかなぁ？
"縦n個、横n個のマス目のそれぞれに 1,2,3,...,n の数字を入れていく。
このマス目の横の並びを行といい、縦の並びを列という。
どの行にも、どの列にも、2つの対角線上にも同じ数字が1回しか現れない入れ方は何通りあるか求めよ。(2020京大文系 改)
n=4がオリジナル
"
library(gtools)
n=4
(a=permutations(n,n)) # ｎの順列
r=nrow(a)
f<-function(x){ # x=c(1,2,3,4) -> rbind(a[1],a[2],a[3],a[4])
ans=NULL
for(i in x){
ans=rbind(ans,a[i,])
}
return(ans)
}
b=permutations(r,n)
head(b) ; tail(b)
B=apply(b,1,f)
diag2 <- function(x){ # 行列xの右上から左下への対角線の配列を返す
n=nrow(x)
ans=numeric(n)
for(i in 1:n) ans[i]=x[i,n+1-i]
return(ans)
}
g <- function(v){ # matri(v,n)で列、対角線の要素がすべて異なればTRUEを返す
n=sqrt(length(v))
(x=matrix(v,n))
if(length(unique(diag (x))) < n) return(FALSE)
if(length(unique(diag2(x))) < n) return(FALSE)
flag=TRUE
for(i in 1:n){
if(length(unique(x[,i])) < n){
flag=FALSE
break
}
}
return(flag)
}
counter=NULL
for(i in 1:ncol(B)){
if(g(B[,i])) counter=append(counter,i)
}
length(counter)
counter
matrix(B[,counter[1]],n)
matrix(B[,counter[2]],n)
matrix(B[,counter[48]],n)
n=5 にするとメモリー不足でエラーがでた。
サンプリングでやってみる
# simulation
n=5
sim <- function(n=5){
v=as.vector(t(replicate(n,sample(n)))) # n × n行列要素をベクトル化
if(g(v,dia=TRUE)){ # diagonal latin squareなら
print(matrix(v,n)) # 行列化して表示
return(1)
}
return(0)
}
k=1e6
sum(replicate(k,sim(5)))
[,1] [,2] [,3] [,4] [,5]
[1,] 2 5 4 3 1
[2,] 4 3 1 2 5
[3,] 1 2 5 4 3
[4,] 5 4 3 1 2
[5,] 3 1 2 5 4
https://id.fnshr.info/2019/07/20/r3-6-0/
で知った
例
rep("1:6",5)
paste(rep("1:6",5),collapse=',')
str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))
eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')))
(s=rep("1:6",5))
(str=paste(s,collapse=','))
(lang=str2lang(paste("expand.grid(",str,')')))
eval(lang)
https://journals.lww.com/investigativeradiology/Abstract/publishahead/The_Clinical_and_Chest_CT_Features_Associated_with.98832.aspx
の付表にこんなデータがあった
# Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value
# Number of lobes Median (interquartile range, IQR)
# 5 (4-5) 5 (5-5) 5 (3-5) 0.003
重症群２５例 中等症群５８例でCTで肺炎像が確認できた肺葉数(１～５）が
全体で中央値が５ 第一四分位数(少ない方から数えて25%の値）が４、第三四分位数(少ない方から数えて75%の値）が５
重症群では中央値５、第一四分位数５、第三四分位数５
中等症群では中央値５、第一四分位数３、第三四分位数５で
p値が０．００３であるという。
０．００２５～０．００３５までとして、
この条件に見合う８３例をリストアップすると
> re[100,][1:25]
[1] 1 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
> re[100,][26:83]
[1] 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[48] 5 5 5 5 5 5 5 5 5 5 5
という類の組み合わせで
19504通りと計算された。
重症群の肺葉数の平均は４．７６
中等症群の平均は４．０４
と計算された。
タイ値がこれだけ多いデータにWilcox検定を使う意味がわからんな。
# アルゴリズムには自信がないｗ
sim <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){
# n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画
prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定
PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算
prev.post = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算
# NPVでも同じ結果
# NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算
# prev.post = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算
mean=mean(prev.post) # 期待値
median=median(prev.post) # 中央値
mode=density(prev.post)$x[which.max(density(prev.post)$y)] # 最頻値
LU=unlist(HDInterval::hdi(prev.post))[1:2] # 95%信頼区間
re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値
if(print){ # ヒストグラム描画
par(mfrow=c(2,2))
hist(prev,freq=F,xlim=c(0,1)) ; lines(density(prev),lwd=2)
hist(PPV,freq=F,xlim=c(0,1)) ; lines(density(PPV),lwd=2)
hist(prev.post,freq=F,xlim=c(0,1)) ; lines(density(prev.post),lwd=2)
BEST::plotPost(prev.post)
par(mfrow=c(1,1))
}
print(round(re,4)) # 概算値表示
invisible(list(re,prev.post)) # 代表値と事後有病率を非表示で返す
}
sim(100,1)
sim(1000,10)
sim(10000,100)
sim(100000,1000)
間違いに気づいたので撤回します。
【富山県最強伝説】新型コロナウイルスPCR検査件数 54人 陽性0人
ある集団から54人を無作為に選んでPCR検査したら陽性0であった。
PCR検査の感度0.7 特異度0.9としてこの集団の有病率の期待値と95%信頼区間を推測せよ。
#----------------------------------------------------------------------------
# 感度SEN, 特異度SPCの検査でN人中X人が陽性であったときの推定有病率prevalence
# 弱情報事前分布はprevalence ~ dunif(0,UL) , UL:上限
#----------------------------------------------------------------------------
PCRj <- function(N,X,UL=1,SEN=0.7,SPC=0.9,verbose=TRUE){ # UL:upper limit of dunif(0,UL)
library(rjags)
modelstring=paste0('
model
{
x ~ dbin(p,n)
p <- prev*sen + (1-prev)*(1-spc)
prev ~ dunif(0,',UL,')
}
')
if(verbose & UL!=1) cat(modelstring)
writeLines(modelstring,'TEMPmodel.txt')
dataList=list(n=N,x=X,sen=SEN,spc=SPC)
jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList,quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p"), n.iter=1e6, thin=10)
js=as.matrix(codaSamples)
BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) ; lines(density(js[,'prev']),col='skyblue')
round(c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]),10)
}
PCRj(100,1)
PCRj(1000,10)
PCRj(10000,100)
PCRj(100000,1000)
PCRj(54,0)
PCRj(54,0,UL=0.1)
1 -0.0377710 0.02994835 0.1779426 0.29057472 0.13252796 0.3279709 -0.3534731 -0.54812874 -0.3869768 -0.1233967
2 0.2707033 0.09238235 -0.1042958 -0.08485129 -0.02590955 0.2159095 -0.3338677 0.00865952 -0.1575856 0.0953795
library(psych)
ICC(DF）
で級内相間求めようとすると変な値でるんだけどなんででしょう？
V1 V2
a -0.03777100 0.27070330
b 0.02994835 0.09238235
c 0.17794260 -0.10429580
d 0.29057472 -0.08485129
e 0.13252796 -0.02590955
f 0.32797090 0.21590950
g -0.35347310 -0.33386770
h -0.54812874 0.00865952
i -0.38697680 -0.15758560
j -0.12339670 0.09537950
> library(psych)
> ICC(x=DF,lmer=FALSE)
Call: ICC(x = DF, lmer = FALSE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.36 2.1 9 10 0.13 -0.18 0.74
Single_random_raters ICC2 0.34 2.0 9 9 0.17 -0.26 0.74
Single_fixed_raters ICC3 0.32 2.0 9 9 0.17 -0.24 0.72
Average_raters_absolute ICC1k 0.53 2.1 9 10 0.13 -0.43 0.85
Average_random_raters ICC2k 0.51 2.0 9 9 0.17 -0.69 0.85
Average_fixed_raters ICC3k 0.49 2.0 9 9 0.17 -0.63 0.84
Number of subjects = 10 Number of Judges = 2
すいません
書き方が良くなかったです。
subject が１，２
judge a,b,c,dをきちんとfirst,second,third,fourth
と書くべきでした。
subjectとjudgeはあってると思うのです。
first second third fourth
0 0.05 0.1 0.15
0 0 0.05 0.05
0.15 0.1 0.05 0.1
0 0.05 0.2 0.05
0.1 0.15 0.1 0.05
0.15 0.05 0.1 0.05
というデータで
Call: ICC(x = over_sight)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.0068 1 5 18 0.43 -0.19 0.48
Single_random_raters ICC2 0.0068 1 5 15 0.44 -0.19 0.48
Single_fixed_raters ICC3 0.0068 1 5 15 0.44 -0.19 0.48
Average_raters_absolute ICC1k 0.0268 1 5 18 0.43 -1.70 0.79
Average_random_raters ICC2k 0.0268 1 5 15 0.44 -1.82 0.79
Average_fixed_raters ICC3k 0.0268 1 5 15 0.44 -1.82 0.79
Number of subjects = 6 Number of Judges = 4
結果のとり得る値が1～0なので最大15％ぐらいの変動幅で、
制度としてはかなりいい検査方法だとおもったのですが・・・
範囲制約性の問題があるのですね。
今回のサンプルがみんな低い得点で分散が低いからICCが低く算出されていました。
ダミーデータ
1,1,1,1
を追加したらICC=0.99となりました。
得点の高いサンプルを追加しないとICCで信頼性は測れないのですね。
上記のような場合で、健常者のみでも信頼性を算出する方法っていうのはないのでしょうか？
DF=rbind(
c(0,0.05,0.1,0.15),
c(0,0,0.05,0.05),
c(0.15,0.1,0.05,0.1),
c(0,0.05,0.2,0.05),
c(0.1,0.15,0.1,0.05),
c(0.15,0.05,0.1,0.05),
c(0.0001,0.0001,0.0001,0.0001))
とダミーに0.0001を加えても
> DF2ICC1(DF)
ICC(1,1) ICC(1,k)
0.2462707 0.5665262
> psych::ICC(DF)
type ICC F df1 df2 p lower bound
Single_raters_absolute ICC1 0.25 2.3 6 21 0.072 -0.027
Single_random_raters ICC2 0.25 2.3 6 18 0.079 -0.028
Single_fixed_raters ICC3 0.25 2.3 6 18 0.079 -0.034
Average_raters_absolute ICC1k 0.57 2.3 6 21 0.072 -0.115
Average_random_raters ICC2k 0.57 2.3 6 18 0.079 -0.122
Average_fixed_raters ICC3k 0.57 2.3 6 18 0.079 -0.154
となるから、高値データでなくてもいいみたい。
https://cran.r-project.org/doc/manuals/r-devel/NEWS.html
install.packages(available.packages()[,1])
これやってからでないとパッケージのコンパイルができない。
After installation is complete, you need to perform one more step to be able to compile R packages:
you need to put the location of the Rtools make utilities (bash, make, etc) on the PATH.
The easiest way to do so is create a text file .Renviron in your Documents folder which contains the following line:
PATH="${RTOOLS40_HOME}\usr\bin;${PATH}"
You can do this with a text editor, or you can even do it from R like so:
writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")
libraryのフォルダで dir /B とやって既にinstall済のリストを作って
install.packages("rstan","ggplot2",....)とやった方がいいな。
5000ほどパッケージがあるらしいけど、全部インストールしたらどれくらいの容量が必要なんだろう？
targetPackages <- available.packages()[,"Package"]
newPackages <- targetPackages[!(targetPackages %in% installed.packages()[,"Package"])]
install.packages(newPackages, repos = "http://cran.r-project.org";)
ステーキハウスがテイクアウトを始めたからインセンティブはその足しにしようっと。
コンパイルなしでやってみた。
サイズ：19.6 GB (21,082,530,395 バイト)
ディスク上のサイズ：20.3 GB (21,805,858,816 バイト)
ファイル数: 527,982、フォルダー数: 118,521
自己解決しました。ライブラリの不足でした。
col=番号指定での色が目に優しくなったみたい。
hist(runif(100),col=2)
hist(runif(100),col=rgb(1,0,0,1))
hist(runif(100),col='red')
colours() と入力すると 657色から選べるから256以上ありそう。
The palette() function has a new default set of colours (which are less saturated and have better accessibility properties).
There are also some new built-in palettes, which are listed by the new palette.pals() function.
These include the old default palette under the name "R3". Finally, the new palette.colors() function allows a subset of colours to be selected from any of the built-in palettes.
https://developer.r-project.org/Blog/public/2019/11/21/a-new-palette-for-r/index.html
binom.ci <- function(x,n,cl=0.95){
PEci=binom::binom.confint(x,n,conf.level=cl)
y=nrow(PEci):1
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mar=c(3,3,2,2))
par(oma=c(0,4,0,0))
plot(PEci$upper,y,xlim=c(min(PEci$lower),max(PEci$upper)),type="n",
yaxt="n",ylab="",xlab="Confidence Interval",
main=paste("C.I. for", x,"out of",n))
points(PEci$upper,y,pch="|") ; points(PEci$lower,y,pch="|")
segments(PEci$lower,y,PEci$upper,y,lwd=3,
col=sample(colours(),1))
abline(v=PEci$mean[1],lty=3)
axis(2, at=y, labels=PEci$method,as.vector,las=2,cex=0.75)
PEci
}
binom.ci(3,312)
情報、ありがとうございます。
?paletteのヘルプにあるスクリプトを実行すると体感できました。
## Demonstrate the colors 1:8 in different palettes using a custom matplot()
sinplot <- function(main=NULL) {
x <- outer(
seq(-pi, pi, length.out = 50),
seq(0, pi, length.out = 8),
function(x, y) sin(x - y)
)
matplot(x, type = "l", lwd = 4, lty = 1, col = 1:8, ylab = "", main=main)
}
sinplot("default palette")
palette("R3"); sinplot("R3")
palette("Okabe-Ito"); sinplot("Okabe-Ito")
palette("Tableau") ; sinplot("Tableau")
palette("default") # reset
theta[i,j] <- min(1, exp(-alpha[i]*t[j]) + beta[i])といいれて、確率を１以下にして二項分布させてようとすると
Node inconsistent with parents というエラーがでた。 stackoverflow.comも明確な答が見いだせず、あれこれいじって
どうも確率1や0で二項分布する可能性があるとエラーがでるみたい、というのを発見。
1を0.99999とするとエラーが回避できた。
1を 1 - 1e-16としてもエラー回避できた。
ある人物Dが新型コロナ肺炎に罹患したとする。
行動調査によって発症前にキャバクラに行っており
接客したキャバ嬢Cが人物D発症の２日後に発症していたことがわかった。
キャバ嬢Cは人物Dから移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢Cから移された可能性もあると主張して
その確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
"
潜伏期間が従う分布はこれ
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612
Gt <- function(delay){
C=rlnorm(1e6,ln_par1,ln_par2)
D=rlnorm(1e6,ln_par1,ln_par2)
mean(C-D > delay)
}
Gt(2)
発生させた乱数から確率密度関数を作ってくれる。
複雑な分布だと収束できず作れないこともあるけど
対数正規分布の差の分布を乱数発生させてpdfが作れた。
library(polspline)
delay=2
ln_par1 = 1.434065
ln_par2 = 0.6612
c=rlnorm(1e5,ln_par1,ln_par2)
d=rlnorm(1e5,ln_par1,ln_par2)
hist(c-d, freq=F,col='white',breaks=100,ylim=c(0,0.11),
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
fit=logspline(c-d)
curve(dlogspline(x,fit),-30,30,ann=F,bty='n',add=T)
segments(delay,0,delay,dlogspline(delay,fit),pch=19,col=4)
curve(dlogspline(x,fit),delay,30,add=T,type='h',col=4)
curve(plogspline(x,fit),-30,30,lty=3,bty='n')
1-plogspline(delay,fit)
fn <- function(delay) 1- plogspline(delay,fit)
curve(fn(x),0,14, bty='n' ,xlab='Delay', ylab='Probability')
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/36.html
mytest <- function() {
ANS <- readline("1+2? ")
if (ANS == "3") cat("Correct!\n")
else cat("Wrong...\n")
}
mytest()
mymessage <- function() {
Sys.sleep(3)
cat("Hello.\n")
}
mymessage()
尚、GIF作製には
https://www.screentogif.com/
がお手軽だった。
'%&%' <- function(x,y) paste0(x,y)
という関数が定義されていた。
これは便利と思った。
> '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi
[1] "√2 =1.4142135623731, π = 3.14159265358979"
でベクトルの内積とかも
'%.%' <- function(x,y) sum(x*y)
%>%パイプとの親和性が無いのが痛い
それって意味があるの？
> paste0('√2 =', sqrt(2), ', π = ', pi)
[1] "√2 =1.4142135623731, π = 3.14159265358979"
普通にpaste0のまま使った方がタイプ量が圧倒的に少ないが。
括弧の対応で混乱しそうなときにちょっと役にたつかも。
こんな感じ
eval(str2lang("expand.grid(" %&% paste0(rep("1:2",6),collapse=',') %&% ")"))
うーん、その例だとイマイチに思える
rerun(6, 1:2) %>% expand.grid
で済むものをわざわざ複雑にしているだけに感じる
無理やり、例示した感は否めないな。
標準装備でなら、
expand.grid(replicate(6,1:2,simplify = FALSE))
ですむし。
特にカンマが文字列に含まれるときにはpaste()より見やすくなるはず
sprintf("√2 = %s, π = %s", sqrt(2), pi)
glue("√2 = {sqrt(2)}, π = {pi}")
確かに便利ですね。
カンマを ','で入力するのは混乱のもとでしたから。
有用な情報、ありがとうございます。'
> glue::glue("√2 = {sqrt(2)}, sin(π/3) = {sin(pi/3)}")
√2 = 1.4142135623731, sin(π/3) = 0.866025403784439
ある関数を作りたいのですが、その関数に入れる数値の数はランダムな場合はどのように書けばいいのでしょうか。
例えば、入れた要素をすべて掛け合わせる関数が作りたい場合、そこに入れる要素が（２，３，４）でも、（２，３，４，５，６）でも反応するような関数の書き方を教えてくださると助かります。
そういう話じゃない？
847さん、回答ありがとうございます。
確かにベクトルで何かするのはわかるのですが、その部分がわからないのです。
出来れば下に問題を書くので、コードを書いていただけると非常に助かります。
問、正の整数であるベクトルを入力すると、そのベクトルの要素をすべて掛け合わせた値を結果として返す関数を作りなさい。
問題文これだけ？他に条件はないの？
標準の関数にあるけどいいのかな？
関数fを作るとして、一番いいのは標準の関数に余計な変更を加えないことだから、これが答えだけど…
f <- prod
実は問２がありまして、たぶんそれは既存の関数にはないと思うので、そちらも教えていただけると非常に助かります。
問２、要素がすべて整数であるベクトルを入力すると、そのベクトルの要素の絶対値をすべて掛け合わせた値を結果として返す関数を作成せよ。
f = function (x) abs(prod(x))
while やfor もしくは再帰関数を使って作らせる問題かな？
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans
}
> f(2,3,4)
[1] 24
> f(2,3,4,5,6)
[1] 720
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans=ifelse(ans>0,ans,-ans)
ans
}
> f(-1,2,3)
[1] 6
> f(1 -2,-3,-4,-5)
[1] 60
入力はベクトルだったから、
> f <- function(v){
+ ans=1
+ for(i in v) ans=ans*i
+ ans=ifelse(ans>0,ans,-ans)
+ ans
+ }
> f(c(-1,2,3))
[1] 6
> f(c(1 -2,-3,-4,-5))
[1] 60
f <- function(...){
v=c(...)
n=length(v)
sub<-function(n){
if(n==0) return(1)
else v[n]*Recall(n-1)
}
sub(n)
}
> f(2,3,4)
[1] 24
> f(2,3,4,5,6)
[1] 720
854, 855, 856を真似るなよ。
rec.prod <- function(x) {
# 再帰により x の要素の絶対値の積を求める。
# 実際には prod(abs(x)) が良い。
l <- length(x) # ベクトルの要素数
if (l == 0) # 空のベクトルなら
return(1) # 積は1
x1.a <- abs(x[1])
if (l == 1)
x1.a
else
x1.a * rec.prod(x[-1])
}
皆様のお力添えのおかげで無事問題を解くことができました。
初心者質問にやさしく回答して頂いて感謝しかありません。
本当にありがとうございました。
この方が>856みたいに関数の中で関数を定義とかなくてカッコいいと思ったのでsystem.timeで時間を比較しようと思って実験を始めたら、
何故か、rec.prodの方がoverflowしやすいなぁ。
頭からかけても尻からかけても同じじゃないかと思うんだが。
> prod(1:100)
[1] 9.332622e+157
> rec.prod <- function(x) {
+ # 再帰により x の要素の絶対値の積を求める。
+ # 実際には prod(abs(x)) が良い。
+ l <- length(x) # ベクトルの要素数
+ if (l == 0) # 空のベクトルなら
+ return(1) # 積は1
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rec.prod(x[-1])
+ }
> rec.prod(1:100)
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
> f <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*Recall(n-1)
+ }
+ sub(n)
+ }
> f(1:100)
[1] 9.332622e+157
> rm(list=ls())
> rp <- function(x) {
+ # 再帰により x の要素の絶対値の積を求める。
+ # 実際には prod(abs(x)) が良い。
+ l <- length(x) # ベクトルの要素数
+ if (l == 0) # 空のベクトルなら
+ return(1) # 積は1
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rp(x[-1])
+ }
> f1 <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*Recall(n-1)
+ }
+ sub(n)
+ }
> f2 <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*sub(n-1)
+ }
+ sub(n)
+ }
> system.time(replicate(1e5,prod(1:10)))
user system elapsed
0.11 0.00 0.11
> system.time(replicate(1e5,rp(1:10)))
user system elapsed
2.64 0.00 2.64
> system.time(replicate(1e5,f1(1:10)))
user system elapsed
2.44 0.01 2.47
> system.time(replicate(1e5,f2(1:10)))
user system elapsed
1.86 0.00 1.87
この問題を再帰で解くのはムダ以外の何ものでもない。
喩えれば、
親分が子分を再帰で酷使しているだけだから親分は消耗しない
という理解でいい？
再帰は、複雑な問題をより小さい問題に分割することがキモだ。
856は元のベクトルはそのまま変わらずに存在していて、ある特定の要素を一つずつ処理しているからループとやってることは変わりない。
859は、一つ短いベクトルを自分自身に処理させているから真の再帰なのだ。
f自身は再帰関数じゃないけど、呼び出しているsubは再帰関数では？
こういうのと同じ。
f <- function(x){ # 数列の和の階乗を返す
m=sum(x)
sub <- function(n){
if(n==1) return(1)
else{
return(n*Recall(n-1))
}
}
fact(m)
}
> f(c(1,2,3))
[1] 720
(ご入力訂正)
f自身は再帰関数じゃないけど、呼び出しているsubは再帰関数では？
こういうのと同じ。
f <- function(x){ # 数列の和の階乗を返す
m=sum(x)
sub <- function(n){
if(n==1) return(1)
else{
return(n*Recall(n-1))
}
}
sub(m)
}
> f(c(1,2,3))
[1] 720
856が再帰っぽくないってのは、859はxを分割・再帰していってるのに対して、856はnを減じていってその都度vのn番目の値を参照してるところかな
それから、859がrec.prod(1:13)であっさりオーバーフローしてしまうのは、integer型だから
あと実行速度が少し遅いのは、xl.aへの代入操作と、if(l ==1){}の余分な条件分岐のせいだね
その辺を修正してやれば、こんな感じかな
# 再帰で問２
f <- function(x){
if(length(x) == 0) return(1)
abs(x[1]) * f(x[-1])
}
勉強になりました。
ついでにwhile版を書いてみました。
f <- function(x){
ans=1
while(length(x)){
ans=ans*abs(x[1])
x=x[-1]
}
ans
}
859が遅い根本の原因は、x[-1] が生成されること。
他の指摘は、細かい高速化の工夫だから初心者に薦めるのはどうかと。
現に869で全く実用的に意味のないコードを書いているし。
題意を素直に表したほうがよいのでは？
絶対値の積は、積の絶対値と等しいから、absを呼ぶのは1回で済むが、何も説明無しにそうするのは私的にはアウトだ。メンテナンスの時に自分や他人が苦しむから。
意味不明なこと言いだしてるが大丈夫か、落ち着け
859は、無駄な条件分岐あるし、すぐオーバーフローするし、コードが冗長だしで、それこそ初心者以外にもNGだぞ
絶対値の積うんぬんに関しては、859と同じでabsを毎回呼んでる、余計な代入操作が無いだけだ
落ち着いてコードを見ろ
速度うんぬんの話はこれまで話題に出てきたように、856と比べての話だ、話についてこれないなら無理に入ってこなくていい
別に高速化を目的としてるわけじゃなく、859から冗長な部分やオーバーフローの原因を取り除いた結果高速化されたってだけだ
だが、速度が遅い原因がxl.aへの代入でなく、x[-1]作成が原因てのはそのとおりだ
文章を書いてる途中で、オーバーフローの原因（代入操作はこっちだ）と勘違いしてしまったようだ
これは訂正に感謝だ
おまえこそもちつけ。
再帰は、再帰の停止条件と、自分自身が行う作業に分離できる。
この場合の停止条件は、要素が1つになったときと考えるのが自然だろう。
数学的には要素が0でもよいし、Rの場合はそれでうまくいくから問題はない。だが、そういった詳細を何も説明せずにこれでよい、と示すのはどうかな？
あと、絶対値を取ったのを代入しているのは、この作業が複雑になった場合の修正の確実さを優先するためで、優秀なバイトコンパイラであれば速度の差はなくなる。
それでも、恥ずかしげもなく、一生懸命マウント取ろうとあがく姿は、
滑稽で哀れになってくるわ
の意味が理解できなくて、実験してみた。
> rec.prod <- function(x) {
+ l <- length(x)
+ if (l == 0)
+ return(1)
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rec.prod(x[-1])
+ }
>
>
> rec.prod(1:13)
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
> rec.prod(as.numeric(1:13))
[1] 6227020800
> rec.prod(c(1,2,3,4,5,6,7,8,9,10,11,12,13))
[1] 6227020800
> rec.prod(c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L))
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
integer型のデータだとオーバーフローしちゃうけど、as.numericでdouble型にするとエラーがでない。
一方、
> f <- function(x){
+ if(length(x) == 0) return(1)
+ abs(x[1]) * f(x[-1])
+ }
> f(1:13)
[1] 6227020800
> f(as.numeric(1:13))
[1] 6227020800
> f(c(1,2,3,4,5,6,7,8,9,10,11,12,13))
[1] 6227020800
> f(c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L))
[1] 6227020800
>
とエラーがでないなぁ。
何故だろう？
とメッセージがでてXXXをインストールすると次は
Package ‘YYY’ was installed before R 4.0.0: please re-install it
とでてきてうんざりしていた。
stackoverflow.comで以下の答をみつけて解決した。
.libPaths()
update.packages(ask=FALSE, checkBuilt=TRUE)
"
問題 ： 「シリツ医 ならば (馬鹿 ならば 裏口 である)」という命題と同値な命題はどれか？
1 : シリツ医 ならば (裏口 ならば 馬鹿 である)
2 : 馬鹿 ならば (シリツ医 ならば 裏口 である)
3 : 馬鹿 ならば (裏口 ならば シリツ医 である)
4 : 裏口 ならば (シリツ医 ならば 馬鹿 である)
5 : 裏口 ならば (馬鹿 ならば シリツ医 である)
"
# PならばQ ≡ (P かつ (Qでない))ではない
'%=>%' = function(P,Q) !(P & !Q)
# premise
A = function(S,B,U) S %=>% (B %=>% U)
# choice
B1 = function(S,B,U) S %=>% (U %=>% B)
B2 = function(S,B,U) B %=>% (S %=>% U)
B3 = function(S,B,U) B %=>% (U %=>% S)
B4 = function(S,B,U) U %=>% (S %=>% B)
B5 = function(S,B,U) U %=>% (B %=>% S)
# premise => choice
C1 = function(S,B,U) A(S,B,U) %=>% B1(S,B,U)
C2 = function(S,B,U) A(S,B,U) %=>% B2(S,B,U)
C3 = function(S,B,U) A(S,B,U) %=>% B3(S,B,U)
C4 = function(S,B,U) A(S,B,U) %=>% B4 (S,B,U)
C5 = function(S,B,U) A(S,B,U) %=>% B5 (S,B,U)
# combination grid of TRUE & FALSE
gr=expand.grid(c(T,F),c(T,F),c(T,F))
# test for premise => choice
all(mapply(C1,gr[,1],gr[,2],gr[,3]))
all(mapply(C2,gr[,1],gr[,2],gr[,3]))
all(mapply(C3,gr[,1],gr[,2],gr[,3]))
all(mapply(C4,gr[,1],gr[,2],gr[,3]))
all(mapply(C5,gr[,1],gr[,2],gr[,3]))
D1 = function(S,B,U) B1(S,B,U) %=>% A(S,B,U)
D2 = function(S,B,U) B2(S,B,U) %=>% A(S,B,U)
D3 = function(S,B,U) B3(S,B,U) %=>% A(S,B,U)
D4 = function(S,B,U) B4(S,B,U) %=>% A(S,B,U)
D5 = function(S,B,U) B5(S,B,U) %=>% A(S,B,U)
# test for choice -> premise
all(mapply(D1,gr[,1],gr[,2],gr[,3]))
all(mapply(D2,gr[,1],gr[,2],gr[,3]))
all(mapply(D3,gr[,1],gr[,2],gr[,3]))
all(mapply(D4,gr[,1],gr[,2],gr[,3]))
all(mapply(D5,gr[,1],gr[,2],gr[,3]))
> #.で文字列を分ける
> unlist(strsplit('3.14','.'))
[1] "" "" "" ""
> unlist(strsplit('3.14','\.')
Error: '\.' is an unrecognized escape in character string starting "'\."
>
> unlist(strsplit('3.14','[.]'))
[1] "3" "14"
> unlist(strsplit('3.14','.',fixed=TRUE))
[1] "3" "14"
人に見てもらいたかったら読みやすく書けよ。アンタのプログラムはいつも読みにくい。
備忘録として書いたのでレスがつくとは思ってなかったが、レスthanx.
初めてみるこんなエラーが返ってきた。
roundoff error is detected in the extrapolation table
困ったときのstackoverflow.comで
https://stackoverflow.com/questions/56384330/integrate-function-returning-roundoff-error-after-working-previously
As for the integration problem, package cubature generally does a better job when integrate
と解決法が記載されていた
integrate(fn, lower, upper)$value
でエラーが返ってきたが、
cubintegrate(fn, lower,upper, method='pcubature')$integral
とするとエラーなしに計算してくれた。
（備忘録）
別スレの問題を「プログラムで面積を数値計算させていたらエラー発生
体験したければ読みにいコードの下記をどうぞ。
r= 0.66017210648907 ; s=2.89088405538017
# 面積計算にintegrate で エラー
sub <- function(c,Area1=1){ # c 円Cの中心のｘ座標
b0=(c^2+r^2-s^2)/(2*c) # B(b0,b1) # 円Cと円Oの交点Bの座標
b1=sqrt(r^2-(c^2+r^2-s^2)^2/(4*c^2))
blue = function(x) sqrt(s^2-(x-c)^2) # 円Cの方程式
BLUE= integrate(blue,b0,c+s)$value*2 # 青の面積
green = function(x) sqrt(r^2-x^2) # 円Oの方程式
GREEN=integrate(green,-r,b0)$value*2 # 緑の面積
GREEN+BLUE-Area1 # GREEN+BLUE=1
}
sub=Vectorize(sub)
x=seq(-r-s,r-s,length=100) ; plot(x,sub(x),pch=19) ; abline(h=0,lty=3)
# 面積計算にcubintegrate でエラーなし
library(cubature)
sub <- function(c,Area1=1){ # c 円Cの中心のｘ座標
b0=(c^2+r^2-s^2)/(2*c) # B(b0,b1) # 円Cと円Oの交点Bの座標
b1=sqrt(r^2-(c^2+r^2-s^2)^2/(4*c^2))
blue = function(x) sqrt(s^2-(x-c)^2) # 円Cの方程式
BLUE= cubintegrate(blue,b0,c+s,method="pcubature")$integral*2 # 青の面積
green = function(x) sqrt(r^2-x^2) # 円Oの方程式
GREEN=cubintegrate(green,-r,b0,method='pcubature')$integral*2 # 緑の面積
GREEN+BLUE-Area1 # GREEN+BLUE=1
}
sub=Vectorize(sub)
x=seq(-r-s,r-s,length=100) ; plot(x,sub(x),pch=19) ; abline(h=0,lty=3)
# 定積分
integral <- function(fn,lwr,upr,...){
cubature::cubintegrate(fn,lwr,upr,method='pcubature',...)$integral
}
library(MASS)にareaという面積計算の関数があるのを思い出したのでやってみた。
library(MASS)
r= 0.66017210648907 ; s=2.89088405538017
# 面積計算にMASS::area
sub <- function(c,Area1=1){ # c 円Cの中心のｘ座標
b0=(c^2+r^2-s^2)/(2*c) # B(b0,b1) # 円Cと円Oの交点Bの座標
b1=sqrt(r^2-(c^2+r^2-s^2)^2/(4*c^2))
blue = function(x) sqrt(s^2-(x-c)^2) # 円Cの方程式
BLUE=area(blue,b0,c+s)*2 # 青の面積
green = function(x) sqrt(r^2-x^2) # 円Oの方程式
GREEN=area(green,-r,b0)*2 # 緑の面積
GREEN+BLUE-Area1 # GREEN+BLUE=1
}
sub=Vectorize(sub)
x=seq(-r-s,r-s,length=100) ; plot(x,sub(x),pch=19) ; abline(h=0,lty=3)
こっちはエラーがでない。
return(anova(...))
}
func(fit1, fit2, fit3)
オブジェクトの引数をうけとって、関数内の関数anovaに渡して
ちゃんと動かすにはどうしたらいいですか？
>773みたいに
c(...)もしくはlist(...)で引き出すのでは動作しない？
> mydata <- data.frame(x = runif(10), y = runif(10), z = runif(10))
> xy <- lm(y ~ x, data = mydata)
> zy <- lm(y ~ z, data = mydata)
> func <- function (...) {
return(anova(...))
}
> func(xy, zy)
Analysis of Variance Table
Model 1: y ~ x
Model 2: y ~ z
Res.Df RSS Df Sum of Sq F Pr(>F)
1 8 0.40535
2 8 0.39312 0 0.012234
ちゃんと動くけど。
何がどう問題なのかを明確にして質問しないと助言を得られないよ。
1万までの素数を求める
> rm(list=ls())
> primes <- function(n) (1:n)[-outer(2:n,2:n)][-1]
> system.time(primes(10000))
user system elapsed
0.60 1.92 2.88
> rm(list=ls())
> library(numbers)
> system.time(Primes(10000))
user system elapsed
0.00 0.00 0.03
いつ、どこによくいるかをクラスタリングしたいのですがRは出来るのでしょうか？
SpaTimeClusパッケージ
1 I 82571
2 B 4909
3 B 230
4 B 231
こんな感じのデータフレームを複数ファイルを読み込んで
ファイル別に横につなげる場合どのようにfor文を書けばいいですか?
ありがとう。これで無事横並べられた
https://gist.github.com/ibombonato/02eea96413dd6a7df5cb38054b2c731c
の計算に、一般化して計算できるように関数を書いてみたのだが、
文字列での入力が必要でどうも美しくない。
# 関数y=f(x)の曲線でx=from からx=toまでの長さを求める
CurveLength <- function(f='sin(x)',from=-pi,to=pi){
str=paste("deriv(~ ",f,",","'x',func=TRUE)")
Df <- eval(str2lang(str))
f1 <- function(x) as.numeric(attributes(Df(x)))
f1=Vectorize(f1)
integrate(function(x) sqrt(1+f1(x)^2), from, to, rel.tol = 1e-12)
}
計算としてはあっていると思う。
> CurveLength()
7.640396 with absolute error < 2.6e-13
> CurveLength('log(1+cos(x))',0,pi/2)
1.762747 with absolute error < 2e-14
もっとエレガントなコードがあればご教示をお願いします。
コードのことならプログラム板のほうで訊いたほうがよさそう
function(x) が \(x)と略記可能になった。
ちょっと練習
f <- \(x) sample(100,x,rep=TRUE)
f(50) |> sort() |> unique() |> length()
> f <- \(x) sample(100,x,rep=TRUE)
> f(50) |> sort() |> unique() |> length()
[1] 40
詳細は
https://cran.r-project.org/bin/windows/base/NEWS.R-4.1.0.html
で
(\(n) sample(10,n, replace = TRUE) |> unique() |> sort()) |> Vectorize() -> f.new
f.new(1:5)
f.old <- Vectorize(function(n){
sort(unique(sample(10,n,replace = TRUE)))})
f.old(1:5)
(\(x) sqrt(1-x^2)) |> Vectorize() |> curve()
> (\(x) sqrt(1-x^2)) |> Vectorize()) |> curve()
Error: unexpected ')' in "(\(x) sqrt(1-x^2)) |> Vectorize())"
(\(x) sqrt(1-x^2)) |> Vectorize() -> fn
curve(fn)
1/4円を描画
> (\(x) sample(1000,x,replace=TRUE) |> sort() |> unique() |> length()) -> f.new
> system.time(replicate(1e5,f.new(5000)))
user system elapsed
28.66 0.01 28.70
> f.old <- function(x) length(unique(sort(sample(1000,x,replace=TRUE))))
> system.time(replicate(1e5,f.old(5000)))
user system elapsed
28.65 0.00 28.67
なので実態としては演算子や関数ではないのでオーバーヘッドが生じないとのこと。
レスありがとうございます。早速、やってみました。
> quote((\(n) sample(10,n, replace = TRUE) |> unique() |> sort()) |> Vectorize())
Vectorize((function(n) sort(unique(sample(10, n, replace = TRUE)))))
今回、搭載された|>は()がないとエラーになるな。
> library(dplyr)
> sample(10,10,replace=TRUE) %>% sort
[1] 1 2 4 5 6 7 7 7 7 8
> sample(10,10,replace=TRUE) %>% sort(decreasing = TRUE)
[1] 9 9 7 6 5 5 3 1 1 1
> sample(10,10,replace=TRUE) |> sort
Error: The pipe operator requires a function call as RHS
> sample(10,10,replace=TRUE) |> sort()
[1] 3 4 5 5 6 6 8 8 8 10
> sample(10,10,replace=TRUE) |> sort(decreasing = TRUE)
[1] 10 10 9 8 7 7 6 6 5 2
オンラインで実行できるサイトだと必ずしも最新版のRに対応していないし、外部ライブラリ非対応のところもあるので
自分でパイプを定義の方がいいかも
'%|%' <- function(x,FUN) FUN(x)
> sample(10,10,replace=TRUE) %|% sort %|% unique
[1] 1 3 5 6 7
hogehoge <- function(x){
cat(deparse(substitute(x)),'=',x,'\n')
}
ch='Swiss'
hogehoge(ch)
> hogehoge(ch)
ch = Swiss
iris[any(iris[,3]==c(6.1,5.1)) ,]
いまいちR言語の仕様がわからん。
適当に推測するに、
iris[iris[, 3] %in% c(6.1, 5.1), ]
雰囲気で伝わるかなと思った。
Rの有利な点は表を容易に分割→編集→結合できる点で、Accessの有利な点はクエリを何段にも重ねることで表を視覚的に扱えることだと思うんだけど、どうかな？
要約すると、細かいとこに手が届くのと、そうではないけど直感的なところ。
でも同じように文字で書くなら、Rの方が便利に思えるな。
Officeが使える環境でR使ってるから、Rだけの環境に移行すればどんな不便が発生するか検討中なんだけど、どう思う？
> SQL(例えばAccess)とRの、データ分析に関する違いって何だと思う？
何を言っているのか分からないが、用意されいている分析手法の豊富さが全く違うだろ。
例えば、多重検定の検定法とかSQLにあるの？ノンパラメトリックの多重検定とかちゃんとできるの？
> SQL(例えばAccess)とRの、データ分析に関する違いって何だと思う？
何を言っているのか分からないが、用意されいている分析手法の豊富さが全く違うだろ。
例えば、多重検定の検定法とかSQLにあるの？ノンパラメトリックの多重検定とかちゃんとできるの？
組み合わせて使うもんだよね。SQLとRは。
ただ、境界線をどこで引くべきか微妙なところで悩む。RだとRStudioのレポートとして残せるから、分析手順を残していくことを考えると極力R使った方がいいのかな。
SQLはデータベースのマネジメントシステムで、
Rはデータ分析のツールじゃん
組み合わせて使う意味はない
絶対に使うな
を実行出来た方いますか？
https://twitter.com/5chan_nel (5ch newer account)
で解決しました。
スレを汚してすみませんせした。
excelは商用だけあって機能的には検定回帰成分分析因子分析何でもできてRに劣るところは無いし、統計以外なら最適化が強い
Rより遅いしVBAは醜いし、365なんていきなり環境ぶっ壊されかねないからRに安住するけど
ぜひともRに取り入れたい機能の一つ
https://navaclass.com/mahalanobis/で
ユークリッド距離のzlimは
contour(x,y,z.E,zlim=c(0,3),nlevels = 5,add=T)
のようにzlim=c(0,3)に設定されているのですが、
マハラノビス距離のzlimは
contour(x,y,z.D,zlim=c(0,1),nlevels = 5,labels="", add=T)
のようにzlim=c(0,1)に設定されています。
確かに、ユークリッド距離のほうは(1, 1)の点が1.0の等高線上にありますし、
マハラノビス距離のほうもAは0.75、Bは2.2という距離に相応しいところに点があります。
ただ、zlim=c(0,3)の意味はきっと0から3までの範囲に限るという意味なんでしょうが、
z.Eとの値も見ても3で区切られているようには見えません。
実験してみましたが、z.Eのほうのzlimは3.5、z.Dのほうのzlimは1.4で等高線が変化するようです。
これでは新たな相関係数を使ったときにどう設定すればよいのか分かりません。
どうか教えて下さい。
xlim, ylimは指定した値でスパッと切るわけではないので（前後に余裕を持たせるようになっている）、
zlimも同じ仕様だと想像するよ。
xlimやylimで指定した特定の値でスパッと切るオプションがあったと思うので、
zlimにも有効でないかどうか、そこから調べて見れば？
ありがとうございます。
今、いろいろ計算した結果、等高線の位置自体は正しいことが判りました。
ただ、問題は等高線の刻みがzlimと相関係数の組み合わせによって変わることです。
具体的には、等高線が0.5刻みになったり1.0刻みになったり0.25刻みになったりします。
にもかかわらずnlevelsは5本のままなので、全体が膨らんだり縮んだり、見た目が全然違ってくるんですよね。
その度に手でzlimの値を1とか3とか5に変えてあげないといけません…これらの数字はまったく理解できません。
なので質問を変えさせていただいて、
相関係数が0以外の時でも等高線上に0.5, 1.0, 1.5, ...などのラベルを強制的に表示する方法は無いでしょうか？
相関係数が0（つまりユークリッド距離）の時は等高線上にそれらの数値が表示されます。
それが相関係数が0.1にでもなった途端に表示されなくなります
（数値表示用に輪が途切れた形でブランクは空いているようなのですが）。
等高線上に0.5, 1.0, 1.5, ...などのラベルが表示されるようになれば、上記の問題も許容できると思います。
余談ですが、
>スパッと切るオプション
調べてみると、xaxs="i"のことのようですね。
"r"だとグラフの端に4%の余白があったのが"i"だと無くなりました。
とりあえず、これで問題ないか確認します。
対象データ
データフレームA内のB列, <chr>型 中身は例として 文字列で"2022/02/24"が入っている。
やりたいこと
この文字列をRのDate型に変換したい
やってみたこと。
tidyverseを用いて
A %>% select(B) %>% as.Date()
これだとエラーが出ました。
エラー文は「'.'からクラス"Date"へ変換はできません」
どうやったらこのB列をRのDate型に変換できますでしょうか。
パイプ処理したいならlubridateパッケージが便利です。チートシート（英語）が公開されているので、使い方は確認してみてください。
パッケージ名のスペル間違ってたらごめん。
変換できない理由は、dplyr::select()関数の返り値がdata.frame型で
as.Date()関数の引数がベクトル型でなければならないからです。
なので、
as.Date(A$B)
とすれば、Date型でベクトルな返り値を得ることができます。tidyverseを
使いたい（というかパイプで処理したい）なら
A %>%
.$B %>% # ベクトルとして取り出す
as.Date()
となります。
返り値をdata.frame型で得たい場合は、dplyr::mutate()関数を使って
A %>%
dplyr::select(B) %>%
dplyr::mutate(B = as.Date(B))
とします。
tidyverseは、基本的にデータフレームを入力としデータフレームを出力と
しますので、列をベクトルとして処理したい場合はmutate()関数を使う必要が
あります。
>934様
ご教示ありがとうございます。
参考にまずは試してみます。
selectじゃなくてpullにすれば動くよ
踏み倒し（＝ハイパーインフレ）の２者選択。
今、日本はお金をどんどん刷って紙幣価値を下げ、
後者の道をまっしぐら。
「税金で借金を返せる難易度ランキング」
１位 日本 ２５７％
２位 スーダン ２１０％
３位 ギリシャ ２０７％
４位 エリトリア １７５％
５位 カーボベルデ １６１％
６位 イタリア １５５％
７位 スリナム １４１％
://twitter.com/fujimaki_takesi/status/1505469787174227973?s=20&t=3LeCpoO4fcZVopj2x8b9Dw
日本から始まる世界的株式市場の大暴落
それが最終的な暴落であることがはっきりするや否や、
マ仆レーヤは出現するでしょう。
マ仆レーヤが公に世界に現れるにつれて、
ＵＦＯがとてつもない数で姿を表すでしょう
https://twitter.com/5chan_nel (5ch newer account)
ベクトル
a <- c(1, 2, 3, 3, 5, 5, 7) #1:7と連番にしたいが、一部値が重複している想定です
に対して、重複を修復するために
df <- data.frame(b=a, c=lag(a, default=a[1]-1))
apply(df, 1, function(x){ifelse(x[1]==x[2], x[1]+1, x[1])})
としましたが、applyを使っているのでかなり遅いです。
質問１：applyを使わない高速化できる手段はないでしょうか。
質問２：ifelseがベクトルに作用すると聞いたので（お恥ずかしいw）
ifelse(df$b == df$a,
としようとしたのですがエラーが出ました。
ifelseは異なる２つのベクトル（行数同じとしても）には作用できないのですね？
何をどうしたいのか、具体的に書いてくれ。
何がしたいのか、さっぱりわからん。
c(1, 2, 3, 3, 5, 5, 7) を 1:7 にしたいだけなら、seq(along=a) でよいが、前と同じ値の場合は、前の値に1を足すというのであれば、rle を使えばなんとかなりそう。
すいませんでした。詳しく書くと、
cベクトルにはcsvから読み込んだ、時間データの秒数が整数値で入っています。
(・csvデータには他に、時間データの年、月、日、時、分、秒、と
全て整数値がそれぞれ別の列で入っています。
・このcsvデータは１秒サンプリングです。
ただそのデータは秒がたまに重複しているのです。
重複している時のルールとして、直前の、１秒前の時間の秒が秒列に入っています。
例で挙げたように
1, 2, 3, 3, 5, 5, 7, 8, 9, 9
この重複を訂正したいのです。
そこで考えたのが 939で挙げたスクリプトだったのですが、
あまりにも遅く(1日分のデータ24*60*60行を訂正するのに）
質問いたしました。
b <- unique(a)
でbは
> b
[1] 1 2 3 5 7
になるよ。
回答ありがとうございます。わかりづらくてすいません。
重複除去ではなく、重複のです。
942に書いたように、
秒データが
1,2,3,4,5,6,7,8
とあるべきところが
1,2,3,3,5,5,7,8
のように4秒のところが直前の３秒に、６秒のところが直前の５秒になっています。
これを修正したいのです。
なら
a <- 1:8
で置換しちゃえばいいんじやないの？
てかまず国語を勉強しよう。
ありがとうございます。
確かに国語力ないですね、今から考えると、、、
必要な前提を書ききれていません。
さらに必要な前提は、データが24*60*60秒分全部なく、１、２秒分抜けているのです。
本当に1秒ごとに確実にデータが来ているなら、seq でも使って新たに値を作ればいいのだが、
きっと、途中でデータが来ない場合もあるのだろう。
例としては、c(1, 2, 2, 2, 9, 9, 11) のようなほうがよいのでは?
この場合、c(1, 2, 3, 4, 9, 10, 11) という答えになればよいのだろう。
r <- rle(a)
res <- lapply(seq(along=r$length), function(i) seq(r$values[i], len=r$length[i]))
do.call(c, res)
これでできると思うが、これで遅いなら、重複のところだけを修正するようにすればよい。
ありがとうございます。
参考に考えてみます。
>>948
返す返す申し訳ないです、、
データの特徴は、
１。重複は２つ。
２。データの抜けは想定できない。
1日に５分ほどのこともあれば(付け加えて言うと、１日のデータをファイルに保存している)
１日に1から３秒ほどの時もあります。(この時は、１日のデータをファイルに保存していない時)
隠さずに具体例を挙げて
FromとToを書けば無駄なやり取りをせずにすむ。
そういうのが「ヘタな抽象化」です。
>>954
ありがとうございます。
データが多いので、下手に抽象化してしまいましたが、
今手元のデータ例では、19:15:26のデータから秒列だけ上げると
下記のように19:15:36から秒列が重複することを始めます。
これが終了するのはデータによりけりで不明です。
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 35, 37, 37, 39, 39, 41, 41, 43, 43, 45
comp.before <- function(x) {
# x の NA を直前の値で補完する.
valid <- !is.na(x) # NA でない(有効な)場所
v <- x[valid] # 有効な値
n <- diff(c(which(valid), length(x)+1)) # それを繰り返す回数
rep(v, n)
}
a <- c(1, 2, 2, 2, 9, 9, 11)
s <- seq(along=a) # 1から始まる等差数列
d <- ifelse(c(FALSE, diff(a) == 0), NA, s - a) # s から a の正しい部分を復元するための差
res <- s - comp.before(d) # 重複部分を修正するには、正しい部分の d を引き継げばよい
寝る前だったので、秒の部分だけでサンプルデータコードにしているので、そこは、
適宜、実際のデータに合わせて修正して。
library(tidyverse)
data.frame(
time = c(26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 35, 37, 37, 39, 39, 41, 41, 43, 43, 45)
) %>%
dplyr::mutate(diff = time - dplyr::lag(time, default = 0)) %>%
dplyr::mutate(new = dplyr::if_else(diff == 0, time + 1, time))
データ重複が2連続だけというのが前提です。3連続以上の場合は、最後の行の
処理を参考に追加してみてください。
ありがとうございます。
なるほど、if_elseをそのように使うのですね。
条件の中で二つのベクトルを使ってダメだったので、諦めてました。
勉強になりました！！
lag()関数のdefault引数の指定があまりよくないので、
こちらの方が良いでしょう。処理手順を分かりやすく
するために処理を分割しています。
ただし、元データに欠損値（NA）がある場合は、意図した
結果を得られない場合がある点には注意してください。
data.frame(
time = c(26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 35, 37, 37, 39, 39, 41, 41, 43, 43, 45)
) %>%
dplyr::mutate(lag = dplyr::lag(time)) %>%
dplyr::mutate(diff = time - lag) %>%
dplyr::mutate(new = dplyr::case_when(diff == 0 ~ time + 1,
diff != 0 ~ time,
is.na(diff) ~ time))
なるほど、さらにありがとうございます。
case_whenやis.na関数を使って、例外処理をなさっているのですね。
非常に勉強になりました。
確かに私があげていたlag関数でのdefault設定だと例外だというのが後で見てわかりづらいですね。
ありがとうございました！！
すいません、前提をいろいろ省略していて、、、
そこは条件文の重ね書きか、modをつかって解決しようとかんがえてました
のであれば、POSIXタイムで処理すれば分またぎも日またぎも対応できるかと。
2022/2/28 23:59:58
2022/2/28 23:59:59
2022/2/28 23:59:59
2022/3/1 00:00:01
2022/3/1 00:00:02
こんな感じのデータでも3行目を「3/1 00:00:00」に修正できるかと。
データが飛んでたりするなら>>957さんの処理を応用すれば良いかと。
そこは正直迷いました、タイムデータを作るのが先か、秒を修正するのが先か。
まずは、秒修正を先にしたのが、Rでのタイムデータの操作と数値データ操作が
どっちが早いかわからなかったので、まずはとっつきやすい数値データを修正して時間データにすることにしたのです。
>>957
礼を述べるのが遅れてすいません。
データが飛んでるところを修正するのに使用いたします。ありがとうございます！！
rstanarm のstan_glm関数を使って簡単な回帰式の練習をしています。
シグマの事前分布として
prior_aux = exponential(0.0008) を設定したいのですが、
Error: '8e-04' is not a supported link for family 'exponential'.（以下略）
というエラーメッセージが出て止まってしまいます。
引数を（autoscale = TRUE) と変えてみても、
Error in exponential(autoscale = TRUE) :
unused argument (autoscale = TRUE)
が出て止まってしまいますので、どうも（）内を認識してくれないようです。
もし対処の方法をご存じの方がおられましたら、ご教示いただければ
ありがたく存じます。
exponentialがrstanarmからじゃなくて他のパッケージから呼び出されてるとか？
試しにrstanarm::exponential(0.0008)で試してみては？
西から出たおひさまが東に沈む
地球の回転に逆らって高速飛行した時の不具合が出てくる。
多分ガウスは勘付いてたっぽい。
966です。早速ご教示いただきまして、どうもありがとうございました。
ご指摘の通りでした。おかげさまで大変勉強になりました。
重ねて厚く御礼申し上げます。
あるcsvでうまくいったスクリプトが、別のスクリプトではうまくいってないのです。
いま、下記のようなtest1.csvがあり、
time
2022/4/4 19:19:10
2022/4/4 19:19:11
下記のスクリプトを通して読み込み、整形しました。
files<-dir(pattern="test1", "./")
tmp <- do.all(bind, lapply(files, function(x) read.csv(x, header=TRUE, stringsAsFactors=FALSE)))
newDate <-as.POSIXct(tmp[,1])
この時、
typeof[,1]は"character"となり、
newDate[1]は、"2022-04-04 19:19:10 JST"となります。
ところが、本番で別の(複数の)csvを読み込ませると、
typeof[,1]は"character"となるのですが、
newDate[1]は、"2022-04-04 JST"となり時刻が消えます。
最初にうまくいったtest1.csvの時刻データが
上記の別のcsvの時刻データをコピーして作ったにもかかわらずです。
何か考えられる原因、解決策はないでしょうか。
a <- c("2022-04-04 23:42:10")
2つ以上 print してみたら？
それが、as.numericしても値がかわってないのです。なので困惑してます。
では、現象が再現する最小のコードとデータを出してください。
きちんと時刻まで入ってないcsvか不正な時刻のcsvが紛れてる
それでした！！
newDate <-as.POSIXct(tmp[,1])
を
library(lubridate)
newDate<-ymd_hms(tmp[,1],tz="Japan")
で変換したところ、ほぼ全てのデータが時刻まで表示されるようになったのですが、
ごく数個、parse errorとなり、NAに変換されてました。
かくにんすると、元データが、"00:00:00"の時だけ、日付になってました。
助かりました。ありがとうございました。
データフレーム内で、数値データに文字列'error'が混入しているとき、
どうやって値を、直前の値に変換したらよろしいでしょうか。
(例)
test <- c(1,2,1,4,'error',8,9,8,'error','error','error',1)
df <- data.frame(c)
となっているデータを、'error'を、その直前の値にしたいのです。
この場合、df$testを(1,2,1,4,4,8,9,8,8,8,8,1)に変換したい。
(**上記データは計測器がはくcsvファイルが元ですが、計測ミスでerrorがまぎれこむのです。)
comp.beforeは適切な名前ではなく、自分はcomplete.by.prev.valueとしている。
縦結合されたセルを処理するときなんかに便利。
NAしか処理してくれないので`error`がNAに置換されている前提のが前提です。
data.frame(
test = c(1,2,1,4,NA,8,9,8,NA,NA,NA,1)
) %>%
tidyr::fill(test)
fill()のリファレンス
https://tidyr.tidyverse.org/reference/fill.html
>>980 >>981 さんが教えてくださったfill関数を使ってみたいと思いましたが、
すいません！、データフレーム内のerrorをNAに変換する方法も教えていただけないでしょうか！
df$test <- as.numeric(df$test)
tidyr::fill(df, test)
ありがとうございました！
なるほど、簡単な手法でさくっとできるのですね、まだまだRになれないと、と痛感しました。みなさん、ありがとうございました！
979に書いてあるやん
csvを読み込ませたいです。
ただし、このcsv、かき特徴があります。
1行目から62行目 いろいろな条件を記載。csvではなく、テキスト。
63行目から ヘッダー含めてcsvファイル形式。
末尾4行 また色々なことを記載。csvではなくテキスト。
最初の62行目までの読み込みをskipさせるのは
read.csv("test.csv", skip=62)
で対処できましたが、末尾4行をスキップさせるのはどうしたら良いでしょうか？？
はい、読み込む時にやりたいです。
後々、フォルダにある数がわからないcsvファイルを一気に読み込んででーたふれーむを作りたいと考えてますので
https://stat.ethz.ch/pipermail/r-help/2008-January/152807.html
最高です！これ！ほんとに助かりました！！！
前からあったっけ？
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html
ちょっと早くなりました。ありがとう。
g = gray, y = yellowとしたいのにg→grayellowとなってしまいます。
正しく置換を行うのはどうすればよいのでしょうか？
>>997の指摘通りで下のようなことをやっているなら、
> x <- factor(sample(c("g", "y"), 20, TRUE))
> x <- gsub("g", "grey", x)
> gsub("y", "yellow", x)
[1] "greyellow" "greyellow" "greyellow" "yellow" "greyellow" "greyellow"
[7] "greyellow" "yellow" "yellow" "greyellow" "yellow" "yellow"
[13] "greyellow" "greyellow" "greyellow" "greyellow" "yellow" "greyellow"
[19] "greyellow" "yellow"
こうするのではなくて、
> gsub("^y", "yellow", x)
[1] "grey" "grey" "grey" "yellow" "grey" "grey" "grey" "yellow"
[9] "yellow" "grey" "yellow" "yellow" "grey" "grey" "grey" "grey"
[17] "yellow" "grey" "grey" "yellow"
こんな風に正規表現で工夫すればよいのでは。
なお、最も簡単な解決策は、green,yerrowの置換順序をyellow,green逆とにすればよい。
s/green,yerrow/grey,yellow/
s/yellow,green/yellow,grey/
>>998
出来ました、ありがとうございます！
新しいスレッドを立ててください。
life time: 1808日 19時間 42分 19秒
運営にご協力お願いいたします。
───────────────────
《プレミアム会員の主な特典》
★ 5ちゃんねる専用ブラウザからの広告除去
★ 5ちゃんねるの過去ログを取得
★ 書き込み規制の緩和
───────────────────
会員登録には個人情報は一切必要ありません。
月300円から匿名でご購入いただけます。
▼ プレミアム会員登録はこちら ▼
https://premium.5ch.net/
▼ 浪人ログインはこちら ▼
https://login.5ch.net/login.php
レス数が1000を超えています。これ以上書き込みはできません。
ニュース
- YouTubeは地上波6局よりも視聴されている スマートTVの利用実態を報告 [PARADISE★]
- タンス預金50兆円が狙われる！4月1日から銀行口座とマイナンバーが紐付けられる『口座管理法』とは…トラブル続出制度に新たな火種 [少考さん★]
- 【芸能】森七菜、サッポロ一番はみそ一択 「母子家庭なので弟に食べさせてあげたら満足してくれた」 [冬月記者★]
- 【NHK】「ママチャリとトラブル」動画を“さらし”拡散 法的責任問われる可能性 ★5 [Hitzeschleier★]
- 【NHK】「ママチャリとトラブル」動画を“さらし”拡散 法的責任問われる可能性 ★4 [Hitzeschleier★]
- 立花孝志氏 若者が選挙に行かないのは「日本が平和だから」自民党政策を称賛「80年も戦死者出してない」 [煮卵★]
- 【WOWOW】UEFAチャンピオンズリーグ・ヨーロッパリーグ 決勝トーナメント★20
- 【WOWOW】UEFAチャンピオンズリーグ・ヨーロッパリーグ 決勝トーナメント★21
- ドジャースvsナショナルズ 7
- 【D専】
- こいせん祝勝会 全レス転載禁止
- 鹿島アントラーズ実況2024☆16
- 【悲報】昨晩の大地震もう風化する [904880432]
- 【速報】大谷、リーグ最多安打の快挙wwwwwwwwwwwwwwww [207516352]
- (´・ω・｀)ショボン好き来て
- みんなを皆んなって書くアホなんなの？
- (´・ω・｀)アニメのパネルが置いてある温泉に行ってきました
- 【悲報】日本人女性さん、500人以上のアメリカ人子供を「実子誘拐」して日本に連れ去っていたwww北朝鮮の拉致事件より酷いのでは・・・ [426633456]