17/09/04 10:49:35.67 xP4OelQr.net
￥
101:￥
17/09/04 12:08:42.71 xP4OelQr.net
￥
102:￥
17/09/04 12:08:59.34 xP4OelQr.net
￥
103:￥
17/09/04 12:09:16.89 xP4OelQr.net
￥
104:￥
17/09/04 12:09:33.74 xP4OelQr.net
￥
105:￥
17/09/04 12:09:50.50 xP4OelQr.net
￥
106:￥
17/09/04 12:10:05.82 xP4OelQr.net
￥
107:￥
17/09/04 12:10:27.40 xP4OelQr.net
￥
108:￥
17/09/04 12:10:59.15 xP4OelQr.net
￥
109:￥
17/09/04 12:11:16.36 xP4OelQr.net
￥
110:１３２人目の素数さん
17/09/05 06:09:36.65 PCfG056b.net
耳栓をしたら世界が変わってワロタ
111:１３２人目の素数さん
17/09/19 11:14:27.03 MBqwivnT.net
Reduce何て便利な関数が標準装備されていたんだな。
?Reduce
Reduce("+", c(1,2,3,4))
Reduce("+", c(1,2,3,4), accumulate = TRUE)
cumsum(c(1,2,3,4))
factorial(4)
Reduce('*',c(1,2,3,4))
Reduce("*", c(1,2,3,4), accumulate = TRUE)
cumprod(c(1,2,3,4))
Reduce("*", c(1,2,3,4), accumulate = TRUE,right=TRUE)
rev(Reduce("*", c(1,2,3,4), accumulate = TRUE,right=TRUE))
cumprod(rev(c(1,2,3,4)))
112:１３２人目の素数さん
17/09/21 19:43:24.15 KStrV4aS.net
96歳の女性をお看取り。肺も心臓も腎臓も機能低下していたが、年齢を考慮して死因は老衰にした。ふと96歳の平均余命はどれくらいかが気になって平成27年簡易生命表（女）を調べてみた。
URLﾘﾝｸ(www.mhlw.go.jp)こういう表をみると自分で計算してみたくなる。諸関数の定義が厚労省の参考資料で解説されているので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))
}
113:１３２人目の素数さん
17/09/21 19:43:57.06 KStrV4aS.net
>>356
実行すると
> LE(f,96)
[1] 3.4
96歳女性の平均余命は3年以上あるらしい。
死因は多臓器不全にする方が正しかったのだろうか？
> LE(m,100)
[1] 2.1
参考資料１ 生命表諸関数の定義
URLﾘﾝｸ(www.mhlw.go.jp)
に寿命中位数というのが定義されている。
半減期にあたる寿命中位数も計算してみた。グラフであたりをつけて内挿してもとめた。
手計算は面倒なので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を参照。
114:１３２人目の素数さん
17/09/21 20:07:38.28 g2Q4cSI1.net
どうでもいいけど、平成27年は完全生命表の年なのに、
なぜわざわざ簡易生命表を使うの？
115:１３２人目の素数さん
17/09/23 10:32:52.85 w2Jlt/XH.net
>>111
5～6年ほど前にやってたな。
生命表と幾何ブラウン運動に基づく株価のモンテカルロシミュレーションで、今の貯蓄で死ぬまでにお金を使い切る確率を計算するの。
116:１３２人目の素数さん
17/09/23 19:46:43.88 lDVqqHRb.net
>>113
公表されていればurl希望
117:１３２人目の素数さん
17/09/23 21:22:11.76 Kb5kBoaF.net
>>115
第22回生命表（完全生命表）
URLﾘﾝｸ(www.mhlw.go.jp)
118:１３２人目の素数さん
17/09/23 21:24:44.71 Kb5kBoaF.net
なお、最新版の平成28年の簡易生命表はこちら
URLﾘﾝｸ(www.mhlw.go.jp)
119:１３２人目の素数さん
17/09/24 06:09:40.57 5AI1PnF/.net
>>117
>>116
ありがとうございました。
120:１３２人目の素数さん
17/09/24 11:13:15.02 o6BLfmLi.net
>>116
## 第22回生命表（完全生命表）
#URLﾘﾝｸ(www.mhlw.go.jp)
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))
}
121:１３２人目の素数さん
17/09/24 11:15:06.76 o6BLfmLi.net
>>117
最新版の平成28年の簡易生命表はこちら
f=c(198,29,19,12,9,7,7,6,5,6,6,7,7,7,8,9,11,12,13,14,16,19,22,24,25,25,25,25,25,26,28,29,31,34,37,41,45,48,50,54,59,66,72,78,84,92,101,112,123,136,149,162,175,186,196,207,220,237,256,275,294,313,334,360,
389,423,463,504,547,595,649,707,776,855,940,1043,1164,1308,1479,1667,1874,2102,2363,2651,2955,3265,3567,3874,4186,4484,4731,4908,5039,5067,4936,4635,4209,3697,3139,2574,2037,1554,1141,805,545,846)
m=c(194,31,21,14,10,9,9,8,7,7,7,7,8,10,13,17,21,26,31,38,44,48,50,51,51,51,53,53,54,56,57,58,60,63,66,70,74,79,83,89,97,106,116,127,140,155,171,190,211,233,255,278,302,328,359,395,434,478,525,572,622,
678,741,813,890,973,1062,1148,1233,1323,1419,1522,1633,1752,1874,2014,2176,2357,2553,2762,2985,3221,3459,3680,3875,4031,4123,4156,4123,4023,3874,3643,3349,3006,2629,2235,1843,1470,1131,837,593,401,258,157,89,90)
LE <-function(ndx,Y,N0=10^5){
n=length(ndx)
lx=numeric(n)
lx[1]=N0
for(i in 1:(n-1))
lx[i+1] <- lx[i] - ndx[i]
nqx=ndx/lx
nLx=numeric(n)
for(i in 1:n)
nLx[i] <- mean(c(lx[i],lx[i+1]))
nLx[n]=0
Tx=rev(cumsum(rev(nLx)))
le=Tx/lx
return(round(le[Y+1],1))
}
LE(f,96)
LE(m,100)
> LE(f,96) # 96歳女性の平均余命
[1] 3.3
> LE(m,100) #100歳男性の平均余命
[1] 1.9
微妙に結果が違うな。
122:１３２人目の素数さん
17/09/27 06:38:20.99 DzhNz5Oe.net
ﾄﾞ底辺特殊シリツ医大に
ｽﾚﾘﾝｸ(hosp板:142番)
のようなスカトロ嗜好の医者がいるかを
1学年100人として無作為の順序で調査していき該当者が1人発見されたら調査終了とするという方法で調査する。
1例目が該当者であった。
この学年のスカトロ嗜好数の期待値を求めよ。
123:１３２人目の素数さん
17/09/27 06:39:04.84 DzhNz5Oe.net
lottery <- function(.N,.r,k=10^4){
f <-function(S,N=.N){
y=c(rep(1,S),rep(0,N-S))
Y=sample(y)
return(Y)
}
g <- function(y){
n=length(y)
for(i in 1:n){
if(!y[i]) i=i+1
else break
}
return(i)
}
xx=0:.N
SS=NULL
for(i in 1:k){
M=t(sapply(xx,f))
for(j in 1:.N){
if(g(M[j,])==.r) SS=c(SS,j-1)
}
}
print(quantile(SS,c(.025,.05,.50,.95,.975)))
c(mean=mean(SS))
}
lottery(100,1)
124:１３２人目の素数さん
17/09/27 06:39:43.99 DzhNz5Oe.net
> lottery(100,1)
2.5% 5% 50% 95% 97.5%
16 22 70 97 98
mean
66.31167
125:１３２人目の素数さん
17/10/02 22:03:04.10 jTVePLTX.net
lottery <- function(.N,.r,k=10^3){
f <-function(S,N=.N){
y=c(rep(1,S),rep(0,N-S))
Y=sample(y)
return(Y)
}
g <- function(y){
n=length(y)
for(i in 1:n){
if(!y[i]) i=i+1
else break
}
return(i)
}
xx=0:.N
SS=NULL
for(i in 1:k){
M=t(sapply(xx,f))
for(j in 1:(.N+1)){
if(g(M[j,])==.r) SS=c(SS,j-1)
}
}
print(quantile(SS,c(.025,.05,.50,.95,.975)))
c(mean=mean(SS))
}
lottery(100,1)
126:１３２人目の素数さん
17/10/06 08:15:11.45 zRAj2AHB.net
Linux版のRでもWindows版みたいに「入力:青字、出力:赤字」みたいなインターフェースにする設定ってある？
127:１３２人目の素数さん
17/10/06 19:21:19.13 nz6qypJh.net
>>125
coloroutパッケージ入れて、色を割り当てるとできるかも。
アーカイブに入っちゃっているけど、使える。
setOutputColors256()のヘルプを参照。
128:１３２人目の素数さん
17/10/07 05:15:49.03 8eYDySnN.net
version 3.4.2になってた。
129:１３２人目の素数さん
17/10/07 09:49:05.80 dws8gMbB.net
順番を表示させる、*1st,*2nd *3rd
> n=21
> paste(n,ifelse(0<n%%10 && n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th'),sep='')
[1] "21st"
130:１３２人目の素数さん
17/10/07 17:31:13.11 O5mMfwME.net
>>128
それだと、11が11stになっておかしいだろ
> library(toOrdinal)
> toOrdinal(21)
[1] "21st"
> sapply(1:20, toOrdinal)
[1] "1st" "2nd" "3rd" "4th" "5th" "6th" "7th" "8th" "9th" "10th"
[11] "11th" "12th" "13th" "14th" "15th" "16th" "17th" "18th" "19th" "20th"
131:１３２人目の素数さん
17/10/07 18:11:16.60 lKxjyQqB.net
11,12,13だけ除外かな
132:１３２人目の素数さん
17/10/07 22:26:53.02 W9vNhGzr.net
>>128
なんでそこで && を使うかなあ。&& の意味わかってる？
ifelse には & を使うことになるのが普通だし、そうでなければ if ～ else のほうがよい。
133:１３２人目の素数さん
17/10/08 10:13:32.31 68bFYZIH.net
>>129
ご指摘とパッケージの御教示ありがとうございました。
134:１３２人目の素数さん
17/10/08 10:15:36.51 68bFYZIH.net
>>131
?&で勉強になりました。
御教示ありがとうございました。
135:１３２人目の素数さん
17/10/08 10:55:49.38 LDCLz3bM.net
>>130
toOrdinal のコードをみたら 11,12,13は別扱いしてあった。
if (tmp %in% c("11", "12", "13"))
tmp.suffix <- "th"
136:１３２人目の素数さん
17/10/08 11:36:14.74 LDCLz3bM.net
legendの中で使うべくone-linerにしたかったのだけど、ここでの指摘に従ってデバッグしたら
paste(n,ifelse(as.character(n%%100)%in%c('11','12','13'),'th',ifelse(0<n%%10 & n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th')),sep='')
こんなに長くなってしまった。
> f<-function(n)paste(n,ifelse(as.character(n%%100)%in%c('11','12','13'),'th',
+ ifelse(0<n%%10 & n%%10<4,c('st','nd','rd')[which(n%%10==1:3)],'th')),sep='')
> prod(sapply(0:10000,f) == sapply(0:10000,toOrdinal::toOrdinal))
[1] 1
一応、パッケージtoOrdinalと0～1万で合致を確認。
いろいろとご助言ありがとうございました。 <(_ _)>ﾍﾟｺﾘ
137:１３２人目の素数さん
17/10/08 11:55:19.96 68bFYZIH.net
パチンコのスレにこういう記載があった。
>バカは何故かパチンコで勝てると思い込んでいる
>バカは本来なら勝てるのに遠隔や不正で負けていると思い込んでいる
>バカは10分の1で10円が当たるクジより1万分の1で8000円が当たるクジの方が儲かると思ってる
>算数すらできないバカが必死に守って支えてきたのがパチンコ業界
これを読んでこんな問題を考えてみた。
宝くじAは1000本中100本が当たりで当たりは 100万円の賞金、
宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、
どちらも売り出し価格は同じなので計100本買うことにする。
Aを何本買うときに賞金の期待値が最大になるか、シミュレーションしてみよ。
138:１３２人目の素数さん
17/10/08 17:27:30.39 d4HM/FAz.net
分散が変わる？
139:１３２人目の素数さん
17/10/08 20:42:24.71 LDCLz3bM.net
シミュレーションすると
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の分布は一様分布になると思ったのだが、
こんな分布になってしまって自分でも混乱している。
URLﾘﾝｸ(i.imgur.com)
お知恵を拝借したい。
140:１３２人目の素数さん
17/10/08 22:12:35.95 COe73DTR.net
よくわからん。rbinomとか使えば？
141:１３２人目の素数さん
17/10/08 22:47:20.08 7IiqlruN.net
>>139
俺もそう思った
142:１３２人目の素数さん
17/10/08 23:23:18.27 7IiqlruN.net
自分なら、こんな感じ
> 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)
同じ結果だよ
143:１３２人目の素数さん
17/10/08 23:29:05.07 7IiqlruN.net
A, Bを少なくとも1本含む100本じゃなかったので、
> b <- function(){which.max(sapply(1:100, a))}
は
> b <- function(){which.max(sapply(0:100, a))}
だな。すまん
> a <- function(i){sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)}
> b <- function(){which.max(sapply(0:100, a))}
> res <- replicate(1000, b())
> hist(res,breaks = 20)
結果は同じ
144:１３２人目の素数さん
17/10/09 00:19:33.56 1+Fkporl.net
>139-142
簡潔なコードありがとうございます。
AもBも賞金の期待値がどちらも10万なので
Aを何本買おうが合計の期待値は同じかなと思いつつ、sample関数を使ってシミュレーションしたのだけれど
一様分布にはならないのでシミュレーションの間違い、
もしくは自分が理解していないsample関数の特性を検出しただけなのか、と思っておりました。
rbinomでのシミュレーション結果とヒストグラムが一致して、一様分布になるという思い込みが間違いだと確信できました。
勉強になりました。ありがとうございました。
145:１３２人目の素数さん
17/10/09 00:27:22.71 1+Fkporl.net
>>142
大したことではないけど
0から始まるので本数とindexを補正して
b <- function(){which.max(sapply(0:100, a))-1}
の方が正確だと思います。
146:１３２人目の素数さん
17/10/09 01:00:19.99 iikzqEdl.net
合計の期待値の標準偏差はAの本数が増えると減るね
147:１３２人目の素数さん
17/10/09 01:05:40.45 1+Fkporl.net
宝くじをA,B併せて100本買ったときの賞金の期待値は
Award <- function(i){ #i:Aを買った本数
sum(dbinom(0:i,i,0.1)*(0:i)*100)+sum(dbinom(0:(100-i),100-i,0.01)*(0:(100-i))*1000)
}
sapply(0:100,Award)
> sapply(0:100,Award)
[1] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
....
[86] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
と買い方に無関係にも思えるんだなぁ。
148:１３２人目の素数さん
17/10/09 02:00:05.31 1+Fkporl.net
>>145
>合計の期待値の標準偏差はAの本数が増えると減る
ことを検証してみた。
a <- function(i) sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)
plot(0:100,sapply(0:100,function(x) sd(replicate(10^3,a(x)))))
URLﾘﾝｸ(i.imgur.com)
149:１３２人目の素数さん
17/10/09 06:15:43.03 1+Fkporl.net
別の話題提供。
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
150:１３２人目の素数さん
17/10/09 15:22:52.36 di3eoSv3.net
どうでもいいんだが、
which.max()は、同値の最大が複数あったとき、最初のインデックスを戻すので、
2番目以降の最大が無視されていることについて、その取り扱いを考慮した方が良いと思うぞ。
> which.max(c(0,1,0,0,0,0,0,0,1))
[1] 2
151:１３２人目の素数さん
17/10/09 17:42:12.61 1+Fkporl.net
>>149
ご指摘ありがとうございます。
最大値を与えるAの購入数が複数ある場合を考えてスクリプトを改変してみました。
a <- function(i){sum(rbinom(i, 1, 0.1) * 100) + sum(rbinom(100-i, 1, 0.01) * 1000)}
which.max2 <- function(x){which(x==max(x))-1}
b <- function(){
b1=NULL
b1=c(b1,which.max2(sapply(0:100, a)))
return(b1)
}
b2=NULL
for(j in 1:1000){
b2=c(b2,b())
}
hist(b2,breaks=20)
152:１３２人目の素数さん
17/10/09 19:43:31.10 1+Fkporl.net
length(b2)
とすると1000を超えているから同値の最大が複数あったことがわかりました。
ありがとうございました。
このスレは勉強になるなぁ。
153:１３２人目の素数さん
17/10/09 22:05:49.50 1+Fkporl.net
どうでもいいことだけど、
hist(b2,breaks=20,col=sample(colors(),2))
とすると毎回配色の変わるツートンカラーのヒストグラムがでて面白い。
154:１３２人目の素数さん
17/10/10 18:01:18.51 S7i/7cX1.net
>>126
125です。レス遅くなったけど、やりたいことができました。ありがとう。
155:１３２人目の素数さん
17/10/11 20:55:10.55 DiiGwH6/.net
array(1:48, dim=c(6, 4, 2))
を，3行2列おきに平均して
,,1
[,1] [,2]
[1,] 5 17
[2,] 8 20
,,2
[,1] [,2]
[1,] 29 41
[2,] 32 44
としたいんですが，効率的な方法はないでしょうか。
実際のデータは，10000x10000x400なんです。
156:１３２人目の素数さん
17/10/11 22:16:13.22 bFq2a2U4.net
>>154
>3行2列
さっぱりわからん。
例えば、結果の[1,1,1]の5は、どれとどれの平均？
157:１３２人目の素数さん
17/10/11 22:25:55.39 HJHHXjdz.net
>>155
6行4列の左上の6つ1,2,3,7,8,9の平均5だと思う。
158:１３２人目の素数さん
17/10/11 22:59:20.81 HJHHXjdz.net
>>154
効率的かどうかはわからんが、やってみた
(d=array(1:48, dim=c(6, 4, 2)))
a=3
b=2
l
159:=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))
160:１３２人目の素数さん
17/10/11 23:00:23.97 UQm0Ph1Q.net
>>155
>>154です。
説明が悪くてすみません。
>>156の通りです。
161:１３２人目の素数さん
17/10/11 23:02:02.80 UQm0Ph1Q.net
>>157
ありがとうございます。
手元にデータがないので明日試してみます。
162:１３２人目の素数さん
17/10/11 23:02:32.81 bFq2a2U4.net
なるほど、じゃあその5,8,17,20はこういうことか。
> a <- array(1:48, dim=c(6, 4, 2))
> b <- expand.grid(1:(dim(a)[1] %/% 3), 1:(dim(a)[2] %/% 2))
> f <- function(i, j, k=1){mean(a[(3*i-2):(3*i),(2*j-1):(2*j),k])}
> mapply(f, b[,1], b[,2])
[1] 5 8 17 20
163:１３２人目の素数さん
17/10/11 23:04:06.69 bFq2a2U4.net
あぁ、すでに >>157 が�
164:嘯曹｢ていたorz
165:１３２人目の素数さん
17/10/11 23:06:59.19 HJHHXjdz.net
> d=array(1:48, dim=c(6, 4, 2))
> a=3
> b=2
> l=6
> m=4
> n=2
> re=NULL
> for(k in 1:n){
+ for(j in 1:(m/b)){
+ for(i in 1:(l/a)){
+ re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),k]))
+ }
+ }
+ }
> array(re,dim=c(m/b,l/a,n))
, , 1
[,1] [,2]
[1,] 5 17
[2,] 8 20
, , 2
[,1] [,2]
[1,] 29 41
[2,] 32 44
データは10000x10000x400なら
l=10000
m=10000
n=400
でよくね？
166:１３２人目の素数さん
17/10/12 07:10:38.11 2sFs7Gsb.net
最初の1行で
> d=array(1:(10000*10000*400))
Error: cannot allocate vector of size 298.0 Gb
エラーがでたｗ
167:１３２人目の素数さん
17/10/12 07:12:09.84 2sFs7Gsb.net
正しいエラーｗはこっち
> d=array(1:(10000*10000*400),dim=c(6, 4, 2))
Error: cannot allocate vector of size 298.0 Gb
168:１３２人目の素数さん
17/10/12 07:21:02.01 2sFs7Gsb.net
そもそも、10000x10000x400だと10000行は3行にならないな。最終行を無視して9999行にするのかな？
169:１３２人目の素数さん
17/10/12 14:48:46.91 sG1WSodN.net
>160を参考に>157のforだらけのスパゲティー・ソースを（ちょっと一般化して）改変しました。
（変数名の重複を避けるため>160の変数名は変更してあります。）
> a=3
> b=2
> c=1
>
> l=6
> m=4
> n=2
>
> A <- array(1:48, dim=c(l,m,n))
> B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c))
> f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])}
> re=mapply(f, B[,1], B[,2],B[,3])
> array(re,dim=c(m%/%b,l%/%a,n%/%c))
, , 1
[,1] [,2]
[1,] 5 17
[2,] 8 20
, , 2
[,1] [,2]
[1,] 29 41
[2,] 32 44
>
上級者のコードを読むのは勉強になります。（深謝）
170:１３２人目の素数さん
17/10/12 16:33:51.03 f/ETd6de.net
先輩方，>>154です。
巨細にわたってご教示・ご議論いただきありがとうございました。
>>166のコードを参考にして，無事処理することができました。
実行速度も申し分なく，効率的に作業が進められそうです。
今回は私には目から鱗な解法で，勉強になりました。
質問して良かったと思っています。。本当にありがとうございました。
10,000行の扱いですが，3行ごとに処理して1行無視してもよいし，
4から7行程度で処理してもよく，フレキシブルです。
今回は4行にして余りなしにしました。
171:１３２人目の素数さん
17/10/12 18:43:54.51 2sFs7Gsb.net
処理速度を比べてみた。
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氏のスキルに改めて感服。
172:１３２人目の素数さん
17/10/14 23:02:55.66 WIgy7a2M.net
# 既約剰余類群 (Z/nZ)*
kjrg <- function(n){
nn=0:(n-1)
f=function(x,y) (x*y)%%n
names(nn)=paste0('n',nn)
z=outer(nn,nn,f)
idx=which(gmp::gcd(1:(n-1),n)==1)+1
return(z[idx,idx])
}
kjrg(10)
kjrg(15)
kjrg(24)
> kjrg(24)
n1 n5 n7 n11 n13 n17 n19 n23
n1 1 5 7 11 13 17 19 23
n5 5 1 11 7 17 13 23 19
n7 7 11 1 5 19 23 13 17
n11 11 7 5 1 23 19 17 13
n13 13 17 19 23 1 5 7 11
n17 17 13 23 19 5 1 11 7
n19 19 23 13 17 7 11 1 5
n23 23 19 17 13 11 7 5 1
173:１３２人目の素数さん
17/10/15 13:21:32.63 GlVyhxtn.net
# 既約剰余類群での加減乗除
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==
174: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
175:１３２人目の素数さん
17/10/15 14:03:02.02 GlVyhxtn.net
名前も一致してidenticalなんだな。
> x=1:10
> y=1:10
> identical(x,y)
[1] TRUE
> names(x)=LETTERS[1:10]
> identical(x,y)
[1] FALSE
> names(y)=LETTERS[1:10]
> identical(x,y)
[1] TRUE
> names(y)=letters[1:10]
> identical(x,y)
[1] FALSE
allの方は内容の全体が一致していれば順番は入れ替わってもTRUEなんだな。
> x=1:10
> y=1:10
> all(x,y)
[1] TRUE
> names(x)=LETTERS[1:10]
> all(x,y)
[1] TRUE
> names(y)=letters[1:10]
> all(x,y)
[1] TRUE
> all(x,sample(y))
[1] TRUE
176:１３２人目の素数さん
17/10/16 01:11:00.56 NsB0m+v3.net
beki <- function(x,n,p){ # x^n%%p
if(n==0) return(1)
if(n==1) return(x%%p)
re=numeric(n)
re[1]=x
for(i in 1:(n-1)) re[i+1] = (x*re[i])%%p
return(re[n])
}
p=1009 # フェルマーの小定理
xx=0:(p-1)
n=p-1
y=sapply(xx,function(x)beki(x,n,p))
summary(y)
177:１３２人目の素数さん
17/10/16 07:46:39.05 NsB0m+v3.net
90歳女性のお看取り
>120のコードから女性の平均余命をグラフにしてみた
Age=85:105
(Life_Expectancy=sapply(Age,function(x) LE(f,x)))
plot(Age,Life_Expectancy,pch=19,col=sample(colors()))
abline(h=5,col='gray',lty=2)
URLﾘﾝｸ(i.imgur.com)
平均余命が5年を切るのは92歳なので、老衰にせずに
主治医指定の病名で診断書を作成。
カルテ記載のしっかりした患者の診療は負担なくこなせて（・∀・）ｲｲ!!
178:１３２人目の素数さん
17/10/16 22:16:44.56 ACageh59.net
outerって使えない場面があるよね？
剰余系での加減乗除での実際
> n=5 # prime number
> nn=1:(n-1)
> tasu <- function(x,y) (x+y)%%n
> hiku <- function(x,y) (x-y)%%n # row - col
> kake <- function(x,y) (x*y)%%n
> g=function(x) nn[which(x==1)]
> .M=outer(nn,nn,kake)
> G=apply(.M,2,g)
> gyaku <- function(x) nn[which(G==(x%%n))]
> waru <- function(x,y) (x*gyaku(y))%%n # row / col
> waru(3,2)
[1] 4
> xx=yy=c(0,nn)
> names(xx)=paste0('x',c(0,nn))
> names(yy)=paste0('y',c(0,nn))
> outer(xx,yy,tasu) # x + y
y0 y1 y2 y3 y4
x0 0 1 2 3 4
x1 1 2 3 4 0
x2 2 3 4 0 1
x3 3 4 0 1 2
x4 4 0 1 2 3
> outer(xx,yy,hiku) # x - y
y0 y1 y2 y3 y4
x0 0 4 3 2 1
x1 1 0 4 3 2
x2 2 1 0 4 3
x3 3 2 1 0 4
x4 4 3 2 1 0
> outer(xx,yy,kake) # x * y
y0 y1 y2 y3 y4
x0 0 0 0 0 0
x1 0 1 2 3 4
x2 0 2 4 1 3
x3 0 3 1 4 2
x4 0 4 3 2 1
> X=Y=nn
> outer(X,Y,waru) # WRONG!!
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] NA NA NA NA
[3,] NA NA NA NA
[4,] NA NA NA NA
> a=expand.grid(X,Y)
> b=matrix(mapply(waru,a[,1],a[,2]),length(X))
> rownames(b)=paste0('x',nn)
> colnames(b)=paste0('y',nn)
> b # x / y
y1 y2 y3 y4
x1 1 3 2 4
x2 2 1 4 3
x3 3 4 1 2
x4 4 2 3 1
179:１３２人目の素数さん
17/10/16 23:28:56.68 FS7yMelZ.net
>>174
そりゃあ、outer に渡す関数がまずいんだろ。
outer に渡す関数は、完全にベクトル対応になっていなければならないんじゃないの？
でなければ、わざわざ outer 使うこともなかろう。
180:１３２人目の素数さん
17/10/17 06:58:05.54 ps5i5Wv8.net
>>175
レスありがとうございます。
自作関数から自作関数を呼び出しているのが原因かと思っていましたが、
ベクトル対応していないとはこういうことでしょうか？
加減乗は引数のどちらをベクトルにしても動作
> kake(1:4,1)
[1] 1 2 3 4
> kake(1,1:4)
[1] 1 2 3 4
> tasu(1:4,1)
[1] 2 3 4 0
> tasu(1,1:4)
[1] 2 3 4 0
> hiku(1:4,1)
[1] 0 1 2 3
> hiku(1,1:4)
[1] 0 4 3 2
除算は第２引数をベクトルにすると誤動作
> waru(1:4,1)
[1] 1 2 3 4
> waru(4,1:4)
[1] 4 1
181:１３２人目の素数さん
2017/10/1
182:7(火) 22:43:30.39 ID:yEir4ohi.net
183:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 07:51:52.36 /ES4CbkI.net
>>177
ありがとうございます。
outerのソースをみたら、こんな記述があったので、X,Yを同時にベクトルで与えても動作する関数でないとダメみたいでした。
FUN <- match.fun(FUN)
Y <- rep(Y, rep.int(length(X), length(Y)))
if (length(X))
X <- rep(X, times = ceiling(length(Y)/length(X)))
FUN(X, Y, ...)
184:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 08:28:19.35 /ES4CbkI.net
# greatest common divisor
gcd <- function(x,y) ifelse(x|y,ifelse(x*y,max(intersect((1:abs(x))[abs(x)%%(1:abs(x))==0],(1:abs(y))[abs(y)%%(1:abs(y))==0])),max(abs(x),abs(y))),NA)
非正整数対応すると長くなるなぁ
> gcd(24,15)
[1] 3
> gcd(16,15)
[1] 1
> gcd(5,0)
[1] 5
> gcd(-24,15)
[1] 3
> gcd(-24,-15)
[1] 3
> gcd(0,0)
[1] NA
185:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 09:17:13.05 /ES4CbkI.net
ユークリッドの互除法をつかうとこんな感じ
> 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
186:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 10:07:52.18 tFigTWxc.net
文系研究者のワイ、困惑
ディープな世界やでほんま
187:１３２人目の素数さん
17/10/25 20:38:52.07 A26rYIkF.net
#.n発r中の狙撃手が.N発狙撃するときの命中数を返す
Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
print(summary(SS))
invisible(SS)
}
188:￥
17/10/26 04:11:35.13 jT09z118.net
￥
189:￥
17/10/26 04:11:58.98 jT09z118.net
￥
190:￥
17/10/26 04:12:18.59 jT09z118.net
￥
191:￥
17/10/26 04:12:35.94 jT09z118.net
￥
192:￥
17/10/26 04:12:53.44 jT09z118.net
￥
193:￥
17/10/26 04:13:10.90 jT09z118.net
￥
194:￥
17/10/26 04:13:28.27 jT09z118.net
￥
195:￥
17/10/26 04:13:46.40 jT09z118.net
￥
196:￥
17/10/26 04:14:03.59 jT09z118.net
￥
197:￥
17/10/26 04:14:22.57 jT09z118.net
￥
198:１３２人目の素数さん
17/10/26 18:17:25.90 tHqsRRzo.net
こういう問題を考えるのは楽しいね。
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
199:１３２人目の素数さん
17/11/01 12:40:01.83 UKV+bz4p.net
ベクトルは再利用されるんだな
> 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
200:１３２人目の素数さん
17/11/03 08:27:39.47 7YwYHhPi.net
> x= list(length=5)
> x
$length
[1] 5
> y=list()
> length(y)=4
> y
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
NULL
> z=list(5)
> z
[[1]]
[1] 5
> w = vector(mode="list" , length= 5)
> w
[[1]]
NULL
[[2]]
NULL
[[3]]
NULL
[[4]]
NULL
[[5]]
NULL
201:１３２人目の素数さん
17/11/09 08:41:04.22 SxG1YX2j.net
Brunner-Munzel検定って比較する両群に重なりがないとｐ値がＮＡになるんだが、
これってバ
202:グなのか？ > 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
203:１３２人目の素数さん
17/11/14 22:02:36.50 Jlt4ZA5T.net
突然ですが質問があります。最近Rソフトを利用し始めた者のです。単刀直入に、
ボロノイ図を制作した場合に分割された面積の平均を出力するには、どうすればよいのでしょうか？
204:１３２人目の素数さん
17/11/14 22:06:32.42 Jlt4ZA5T.net
現状としては住所をグーグルマップの緯度経度に変換し、分割まではできています。
宜しければどなたかお力添えをいただけないでしょうか？
205:￥
17/11/15 00:58:17.28 WZuPK5Ir.net
￥
206:￥
17/11/15 00:58:38.82 WZuPK5Ir.net
￥
207:￥
17/11/15 00:58:56.67 WZuPK5Ir.net
￥
208:￥
17/11/15 00:59:14.81 WZuPK5Ir.net
￥
209:￥
17/11/15 00:59:32.78 WZuPK5Ir.net
￥
210:￥
17/11/15 00:59:50.57 WZuPK5Ir.net
￥
211:￥
17/11/15 01:00:09.89 WZuPK5Ir.net
￥
212:￥
17/11/15 01:00:29.29 WZuPK5Ir.net
￥
213:￥
17/11/15 01:00:47.81 WZuPK5Ir.net
￥
214:￥
17/11/15 01:01:09.56 WZuPK5Ir.net
￥
215:１３２人目の素数さん
17/11/15 17:11:14.46 eZCzfwkb.net
>>198
SpatialPolygonsクラスまでできるといると仮定して、
各ポリゴンの中にareaスロットがあるから、
それを抽出して、その平均を取れば良い。
ボロノイ図がSpatialPolygonsクラスになっていないなら、
まずは変換。
216:￥
17/11/15 20:30:49.23 WZuPK5Ir.net
￥
217:￥
17/11/15 20:31:06.31 WZuPK5Ir.net
￥
218:￥
17/11/15 20:31:21.95 WZuPK5Ir.net
￥
219:￥
17/11/15 20:31:40.54 WZuPK5Ir.net
￥
220:￥
17/11/15 20:31:58.15 WZuPK5Ir.net
￥
221:￥
17/11/15 20:32:16.33 WZuPK5Ir.net
￥
222:￥
17/11/15 20:32:32.65 WZuPK5Ir.net
￥
223:￥
17/11/15 20:32:48.32 WZuPK5Ir.net
￥
224:￥
17/11/15 20:33:05.48 WZuPK5Ir.net
￥
225:￥
17/11/15 20:33:21.41 WZuPK5Ir.net
￥
226:１３２人目の素数さん
17/11/15 21:12:29.48 9zgzmHZc.net
ボロノイ図の各領域の面積の平均値って、単に
描画領域の総面積÷点数 じゃないの？
何か暗黙の前提があるの？
227:￥
17/11/16 06:01:20.12 6ldUKvsQ.net
￥
228:￥
17/11/16 06:01:37.22 6ldUKvsQ.net
￥
229:￥
17/11/16 06:01:56.08 6ldUKvsQ.net
￥
230:￥
17/11/16 06:02:14.27 6ldUKvsQ.net
￥
231:￥
17/11/16 06:02:31.67 6ldUKvsQ.net
￥
232:￥
17/11/16 06:02:48.95 6ldUKvsQ.net
￥
233:￥
17/11/16 06:03:06.96 6ldUKvsQ.net
￥
234:￥
17/11/16 06:03:24.98 6ldUKvsQ.net
￥
235:￥
17/11/16 06:03:42.58 6ldUKvsQ.net
￥
236:￥
17/11/16 06:04:00.50 6ldUKvsQ.net
￥
237:１３２人目の素数さん
17/11/17 20:18:17.65 5lHMDqma.net
Stan って 離散量をparameter 扱いできないので
m ? categorical(θ) としても ｍ の分布を出せない。
JAGSだと
m ～ dcat(θ）
とすると ｍの分布がだせる。
この理解でいいかな？
238:￥
17/11/17 20:41:11.51 Y8c01QBt.net
￥
239:￥
17/11/17 20:41:27.83 Y8c01QBt.net
￥
240:￥
17/11/17 20:41:43.62 Y8c01QBt.net
￥
241:￥
17/11/17 20:41:58.17 Y8c01QBt.net
￥
242:￥
17/11/17 20:42:13.43 Y8c01QBt.net
￥
243:￥
17/11/17 20:42:28.14 Y8c01QBt.net
￥
244:￥
17/11/17 20:42:43.14 Y8c01QBt.net
￥
245:￥
17/11/17 20:42:59.36 Y8c01QBt.net
￥
246:￥
17/11/17 20:43:14.19 Y8c01QBt.net
￥
247:￥
17/11/17 20:43:30.82 Y8c01QBt.net
￥
248:１３２人目の素数さん
17/11/20 17:47:00.70 ptVXRSQM.net
以前、ボロノイ図について質問させていただいた初心者です。
質問後色々と試行錯誤したのですが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
と入力しボロノイ図は出せているのですが、
249:１３２人目の素数さん
17/11/20 17:52:19.12 ptVXRSQM.net
summaryと入力しても
function (object, ...)
UseMethod("summary")
<bytecode: 0x06b5eef0>
<environment: namespace:base>
と表示されます。どの様に入力または再入力すればdir.areaが表示させられるでしょうか？
私自身初めてで、全く分からない状態でRを使用しているので、宜しければ分かり易く教えて頂けると幸いです。
250:￥
17/11/20 21:05:10.85 yjl62pX+.net
￥
251:￥
17/11/20 21:05:28.54 yjl62pX+.net
￥
252:￥
17/11/20 21:05:44.92 yjl62pX+.net
￥
253:￥
17/11/20 21:06:01.26 yjl62pX+.net
￥
254:￥
17/11/20 21:06:18.27 yjl62pX+.net
￥
255:￥
17/11/20 21:06:35.26 yjl62pX+.net
￥
256:￥
17/11/20 21:06:52.55 yjl62pX+.net
￥
257:￥
17/11/20 21:07:09.97 yjl62pX+.net
￥
258:￥
17/11/20 21:07:27.71 yjl62pX+.net
￥
259:￥
17/11/20 21:07:44.86 yjl62pX+.net
￥
260:１３２人目の素数さん
17/11/22 10:59:21.45 1eSvd7ei.net
>>242
検証用のデータがないと誰も追試しないと思う。
>154は模擬データがあったのでやる気になった。
261:￥
17/11/22 12:06:07.41 j3iP9uSb.net
￥
262:￥
17/11/22 12:06:26.73 j3iP9uSb.net
￥
263:￥
17/11/22 12:06:45.80 j3iP9uSb.net
￥
264:￥
17/11/22 12:07:02.42 j3iP9uSb.net
￥
265:￥
17/11/22 12:07:21.82 j3iP9uSb.net
￥
266:￥
17/11/22 12:07:41.50 j3iP9uSb.net
￥
267:￥
17/11/22 12:08:00.58 j3iP9uSb.net
￥
268:￥
17/11/22 12:08:22.66 j3iP9uSb.net
￥
269:￥
17/11/22 12:08:39.11 j3iP9uSb.net
￥
270:￥
17/11/22 12:08:57.12 j3iP9uSb.net
￥
271:１３２人目の素数さん
17/11/24 16:46:01.01 S4G8MO/P.net
関数のオプションの取り得る値を全て表示する方法って知らない？
例えばcor関数のmethodオプションの取り得る値は"pearson", "kendall", "spearman"みたいに。
そういうのは無いのかな？
初めて使う関数のオプションいちいちググるの面倒くさい。
272:￥
17/11/24 18:45:27.78 7RwNGOaZ.net
￥
273:￥
17/11/24 18:45:42.78 7RwNGOaZ.net
￥
274:￥
17/11/24 18:45:58.90 7RwNGOaZ.net
￥
275:￥
17/11/24 18:46:16.24 7RwNGOaZ.net
￥
276:￥
17/11/24 18:46:33.09 7RwNGOaZ.net
￥
277:￥
17/11/24 18:46:49.63 7RwNGOaZ.net
￥
278:￥
17/11/24 18:47:03.93 7RwNGOaZ.net
￥
279:￥
17/11/24 18:47:38.91 7RwNGOaZ.net
￥
280:￥
17/11/24 18:47:57.22 7RwNGOaZ.net
￥
281:￥
17/11/24 18:48:16.46 7RwNGOaZ.net
￥
282:１３２人目の素数さん
17/11/24 19:33:11.50 aKG4T9LS.net
>>265
ヘルプ表示じゃ駄目なのかい？
283:１３２人目の素数さん
17/11/24 20:02:53.10 YIwuPhoP.net
>>265
> args(cor)
function (x, y = NULL, use = "everything", method = c("pearson",
"kendall", "spearman"))
NULL
284:１３２人目の素数さん
17/11/24 22:49:54.91 8tVkv6KQ.net
>>276
前にヘルプで見ようと思ったら書いてなかった覚えが…。マニアックな関数だったんで、そんなもんかな～と諦めてた。やっぱ基本はヘルプか。
>>277
argsは見れる関数と見れない関数があるのであんまり役に立たない印象。
285:￥
17/11/25 00:22:59.73 kb2pRtxG.net
￥
286:￥
17/11/25 00:23:21.36 kb2pRtxG.net
￥
287:￥
17/11/25 00:23:46.66 kb2pRtxG.net
￥
288:￥
17/11/25 00:24:03.98 kb2pRtxG.net
￥
289:￥
17/11/25 00:24:22.25 kb2pRtxG.net
￥
290:￥
17/11/25 00:24:46.55 kb2pRtxG.net
￥
291:￥
17/11/25 00:25:07.98 kb2pRtxG.net
￥
292:￥
17/11/25 00:25:29.17 kb2pRtxG.net
￥
293:￥
17/11/25 00:25:50.80 kb2pRtxG.net
￥
294:￥
17/11/25 00:26:12.77 kb2pRtxG.net
￥
295:１３２人目の素数さん
17/11/25 01:57:31.82 sNF3t9bW.net
>>278
>argsは見れる関数と見れない関数があるのであんまり役に立たない
だったら、例が悪いよ。見られない関数の例をあげないと
> ヘルプで見ようと思ったら書いてなかった
こっちだってわざわざ下請け関数のオプションまでは書かないよ。
ヘルプを見て未解決ならソースを見てくれ。
296:￥
17/11/25 06:17:05.29 kb2pRtxG.net
￥
297:￥
17/11/25 06:17:23.25 kb2pRtxG.net
￥
298:￥
17/11/25 06:17:39.26 kb2pRtxG.net
￥
299:￥
17/11/25 06:17:55.24 kb2pRtxG.net
￥
300:￥
17/11/25 06:18:11.86 kb2pRtxG.net
￥
301:￥
17/11/25 06:18:28.49 kb2pRtxG.net
￥
302:￥
17/11/25 06:18:44.12 kb2pRtxG.net
￥
303:￥
17/11/25 06:19:01.27 kb2pRtxG.net
￥
304:￥
17/11/25 06:19:19.56 kb2pRtxG.net
￥
305:￥
17/11/25 06:19:36.14 kb2pRtxG.net
￥
306:１３２人目の素数さん
17/11/26 17:41:27.37 ahD9OIUm.net
rm(list=ls())
で R のメモリーをクリア
307:１３２人目の素数さん
17/11/28 13:00:14.06 kbsEO2Dt.net
先週長々と質問させていただいたものです。その後、試行錯誤して行った結果十解決し
308:ました。 以前は「summary」としか入れていなかったのですが、「dd$summary」と入力したところ 全体のデータが出力されました。それ以後「dd$summary$dir.area」等で出せるようになりました。 大変初期的な場面で躓いておりましたが、現在何とか進んでおります。 色々とお答えくださった方々、とても参考になりました。本当にありがとうございました。
309:１３２人目の素数さん
17/12/05 15:24:41.67 2L5k1k+Np
semを使用中、specifyModel()でパスを入力後に
NOTE: it is generally simpler to use specifyEquations() or cfa()
see ?specifyEquations
NOTE: adding 3 variances to the model
という文章が出てそのままsem.X <- sem(modelA.X, X, N=n)
と入力すると
solve.default(C[ind, ind]) でエラー:
Lapack routine dgesv: システムは正確に特異です: U[2,2] = 0
追加情報: 警告メッセージ:
1: sem.default(ram, S = S, N = N, raw = raw, data = data, pattern.number = pattern.number, で:
S is not positive-definite: expect problems
2: sem.default(ram, S = S, N = N, raw = raw, data = data, pattern.number = pattern.number, で:
The following variables have no variance or error-variance parameter (double-headed arrow):
x1, x2, x3, x4, x5, x6
The model is almost surely misspecified; check also for missing covariances.
と出るのですがどうすれば警告等が出ずに分析が行えるのでしょうか？
310:１３２人目の素数さん
17/12/10 13:11:49.22 cG4z9akW.net
N(=100)回コインをなげてｎ(=5回)以上続けて表がでる確率。
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
mean(replicate(10^5,seqn()))
> mean(replicate(10^5,seqn()))
[1] 0.81085
案外、高い確率になった。
311:１３２人目の素数さん
17/12/10 22:28:22.98 WuWuRvn7.net
ポアソンでやるやつだっけ
312:１３２人目の素数さん
17/12/11 07:47:07.43 EURogHTh.net
pooledVariance <- function(...) {
args = list(...)
n.args=length(args)
ss2=0
df=0
for(i in 1:n.args){
ss2 = ss2 + var(args[[i]])*(length(args[[i]])-1)
df = df + (length(args[[i]])-1)
}
ss2/df
}
effectsize <- function(y1,y2){
diff=mean(y1)-mean(y2)
var=(var(x1)*(length(x1)-1)+ var(x2)*(length(x2)-1))/(length(c(y1,y2))-2)
sd=sqrt(var)
diff/sd
}
library(effsize)
cohen.d()
313:１３２人目の素数さん
17/12/11 07:50:22.91 EURogHTh.net
>>304
単なる二項分布。
コインが５回続けて表がでたら、0.5^5 <0.05なのでイカサマコインといわれちゃいそうなんだが、
１００回やってみると案外、５回表が続くので確率を計算しようと思ったが、解析的にできる頭がないので
シミュレーションしてみた。
314:１３２人目の素数さん
17/12/11 07:55:42.42 EURogHTh.net
>303は
> 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)
315:%in% rbinom(100,1,0.5)は １個ずつ評価されてTRUEが５個返ってくるだけ。 文字列にしてgrepを使うとなんとかなりそうな気がしないでもないのだけど....
316:１３２人目の素数さん
17/12/22 17:21:16.51 U1Bq8bdb.net
すまん、教えてほしいだけど
分析するために初めてＲをインストールしようと思って、このスレのあるように公式サイト行ったら、esetが「JS/Redirector.NAV トロイの木馬」を検知したんだが…；
どうしたらいいだ…
317:１３２人目の素数さん
17/12/22 20:09:03.17 ol6XbcBE.net
そっちじゃなくてCRAN行け
318:１３２人目の素数さん
17/12/22 20:42:31.02 U1Bq8bdb.net
ホームページダメとかわけわかんねぇ・・・とんでもねぇな
ありがとう
319:１３２人目の素数さん
17/12/22 20:50:20.52 uzZ0hejh.net
>>308
エロサイトにアクセスしてないｗ？
320:１３２人目の素数さん
17/12/22 21:27:01.31 RNoKRtjI.net
>>310
CRANはここ
URLﾘﾝｸ(cran.r-project.org)
プロジェクトの方はESETの誤検知っぽいんだよな
321:１３２人目の素数さん
17/12/23 09:28:44.66 n0SBd+bp.net
>>311
してねぇよｗｗ
totalvirusで調べると一件引っかかるし、Redirector検知だから、しょうがないね
322:１３２人目の素数さん
17/12/24 20:42:51.55 CT/NKMd7.net
col=rgb(runif(1),runif(1),runif(1),runif(1))で色指定すると
走らせるたびに色がちがっておもしろい。
hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
323:１３２人目の素数さん
17/12/31 17:23:19.35 14tdpK/Y.net
stanやJAGSのコードでgamma関数を使おうとして
y = gamma(x)
と、書いたらエラーになった。
stanだと y=tgamma(x)、JAGSだとy=exp(loggam(x))で動作した。
324:１３２人目の素数さん
18/01/02 08:46:50.95 qdmBZ37O.net
ある大学の入学者男女の比率は１であるという帰無仮説を検定する課題が花子と太郎に課された。
花子は５０人を調査できたら終了として入学者を５０人をみつけて１８人が女子であるという結果を得た。
帰無仮説のもとで
５０人中１８人が女子である確率は 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
で帰無仮説は棄却できないと結論した。
URLﾘﾝｸ(i.imgur.com)
一方、本番と十八番が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は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
URLﾘﾝｸ(i.imgur.com)
で帰無仮説は棄却される。
どちらの検定が正しいか、どちらも正しくないか？
検定する意図によってｐ値が変わるのは頻度主義統計の欠陥といえるか？
花子の横軸は裏口入学者数、太郎の横軸はサンプル数なので
サンプルでの裏口入学率を横軸にして９５％信頼区間を示す。
花子の検定での信頼区間は0.36～0.72で18/50を含む、p=0.06491
URLﾘﾝｸ(i.imgur.com)
太郎の検定での信頼区間は0.375～0.72で18/50を含まない、p= 0.0275
URLﾘﾝｸ(i.imgur.com)
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。
325:１３２人目の素数さん
18/01/03 10:35:48.54 YJfyxrv+.net
（訂正）
ある大学の入学者男女の比率は１であるという帰無仮説を検定する課題が花子と太郎に課された。
花子は５０人を調査できたら終了として入学者を５０人をみつけて１８人が女子であるという結果を得た。
帰無仮説のもとで
５０人中１８人が女子である確率は 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
で帰無仮説は棄却できないと結論した。
URLﾘﾝｸ(i.imgur.com)
一方、十八という数字が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は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
URLﾘﾝｸ(i.imgur.com)
で帰無仮説は棄却される。
どちらの検定が正しいか、どちらも正しくないか？
検定する意図によってｐ値が変わるのは頻度主義統計の欠陥といえるか？
花子の横軸は女子数、太郎の横軸はサンプル数なので
サンプルでの女子の割合を横軸にして９５％信頼区間を示す。
花子の検定での信頼区間は0.36～0.72で18/50を含む、p=0.06491
URLﾘﾝｸ(i.imgur.com)
太郎の検定での信頼区間は0.375～0.72で18/50を含まない、p= 0.0275
URLﾘﾝｸ(i.imgur.com)
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。
326:１３２人目の素数さん
18/01/15 16:30:41.76 wJofbCL/.net
kainokousiki<-function(a,b,c){return (-b+sqrt(b^2-4*a*c))/(2*a)} #解の公式
kainokousiki(1,-5,6)
でrunすると3じゃなくて6を返すんだけど、どこが間違ってる？
327:１３２人目の素数さん
18/01/15 16:35:00.20 wJofbCL/.net
自己解決
かっこが足りなかった
328:１３２人目の素数さん
18/01/19 12:12:23.00 k3h1PrrK.net
Pythonのスレはないのか
329:￥
18/01/21 20:11:34.35 oUqQkvBY.net
￥
330:￥
18/01/21 20:11:56.34 oUqQkvBY.net
￥
331:￥
18/01/21 20:12:22.56 oUqQkvBY.net
￥
332:￥
18/01/21 20:12:47.97 oUqQkvBY.net
￥
333:￥
18/01/21 20:13:05.14 oUqQkvBY.net
￥
334:￥
18/01/21 20:13:34.57 oUqQkvBY.net
￥
335:￥
18/01/21 20:14:02.66 oUqQkvBY.net
￥
336:￥
18/01/21 20:14:16.99 oUqQkvBY.net
￥
337:￥
18/01/21 20:14:35.47 oUqQkvBY.net
￥
338:￥
18/01/21 20:15:07.01 oUqQkvBY.net
￥
339:１３２人目の素数さん
18/02/05 19:08:34.64 AIJT+apj.net
機械学習をきっかけにPythonに逆転された感じだね
340:１３２人目の素数さん
18/02/06 18:35:39.90 tAZA/Fp/.net
『Rを使った～』だとPythonじゃないのかよって思うよね
341:１３２人目の素数さん
18/02/06 21:37:32.10 tUqX17n9.net
ずっとRだけでPython触ったこと無いけど、覚え直す価値ある？
環境構築からもう面倒なイメージ
342:１３２人目の素数さん
18/02/06 22:43:23.82 qQMNyZjW.net
Python自体は:と直後のインデントさえ気を付ければかなり簡単
3系は数が全て小数扱いなので楽
Anacondaというパッケージでインストールすれば、今流行りのJupyter Notebookという開発環境で対話的にコーディングできる（Rも使える）
343:１３２人目の素数さん
18/02/07 20:40:49.48 mgaw9oVw.net
アナコンダてのがRぽくできるのね、ありがとう
dplyrやggplot2みたいに素人でも簡単便利だといいんだけど
Pythonはオブジェクト志向ぽいしすぐ諦めそう
344:１３２人目の素数さん
18/02/12 16:19:11.21 NSJ4iUa4.net
ブラウザ環境なくなっちゃったの？
345:１３２人目の素数さん
18/02/12 17:06:56.33 hMO/kWPg.net
>>336
誤爆？
何をブラウズする環境が無くなったの？
346:１３２人目の素数さん
18/02/13 18:21:54.45 1hic99Cx.net
ブラウザでプログラミングする環境
347:１３２人目の素数さん
18/02/13 21:09:53.07 EFzybxrC.net
>>338
RStudio Serverの話？
348:１３２人目の素数さん
18/02/13 22:27:26.57 Qyl0hNg7.net
RStudio Severあるで。Dockerで使うのがいいんじゃないかな？！
349:１３２人目の素数さん
18/02/25 18:12:28.68 aD34K55o.net
統計学とウェブ解析を交えて実践的な勉強と練習を
したいのですが、おすすめな書籍やサイトはありますか。
実際に解析ツールや分析ツールを用いて
自分で分析解析してから
解答を見て適切な手順や方法、考察を
解説
350:してくれるものが良いです。 統計学は統計検定2級の知識はありますが ウェブ解析はテキスト読んだだけです。
351:１３２人目の素数さん
18/02/26 12:43:55.13 U8DvzUE2.net
ウェブ解析とは具体的になんですか
352:１３２人目の素数さん
18/02/27 00:54:41.17 hS0OJ3qQ.net
R studioは日本語コメント書く度にIMEが無効になったりカーソルがずれたり黒文字の予測変換が黒背景と重なって見えなくなったりと散々だわ
353:１３２人目の素数さん
18/02/27 01:06:58.09 GkJoHhTZ.net
Windows環境で使うから
354:１３２人目の素数さん
18/02/27 01:16:41.39 hS0OJ3qQ.net
会社規定なんでしかたない
UNIX環境使えるのが羨ましい
355:１３２人目の素数さん
18/02/27 01:25:59.40 34tpClZB.net
Docker for WindowsでRStudioサーバー動かせば？
356:１３２人目の素数さん
18/02/27 12:19:33.31 ZcuB0rsn.net
なぜPythonスレッドがない？
357:１３２人目の素数さん
18/02/27 14:05:09.17 O+8uJ5V+.net
板名読める？
358:１３２人目の素数さん
18/03/04 20:14:10.91 nWbPE6nD.net
>>343
IME無効になるのは俺だけじゃなかったと知ってほっとした。
359:１３２人目の素数さん
18/03/04 21:28:34.58 R7ZPBSuG.net
ストアアプリも同じ症状でるからRStudio固有の問題でなくWindows環境の不治の病だと思ってる
360:１３２人目の素数さん
18/03/26 20:57:59.83 WAyucZm1.net
Run any R code you like. There are over three thousand R packages preloaded.
URLﾘﾝｸ(rdrr.io)
361:１３２人目の素数さん
18/03/26 21:11:25.16 WAyucZm1.net
>>351
日本語あると動作しないが、
グラフを描いてくれるのはうれしい。
URLﾘﾝｸ(imagizer.imageshack.com)
362:１３２人目の素数さん
18/04/04 07:52:09.98 PZp+DZN4.net
Fisher test検定時に
p<2.2e-16
と表示されるんですが、これより小さい値の指数桁数を正確に表記する方法教えて下さい。
例えば5.8e-35となるようにです。
363:１３２人目の素数さん
18/04/04 08:05:49.72 lBVECsOD.net
返値の中にあるp値を参照しなされ
364:１３２人目の素数さん
18/04/04 10:02:52.31 CRvlhZKw.net
fisher.test関数の返り値はlist型で、その中にp.valueという名前でp値が格納されているから$演算子を使って直接参照するか、broom::tidy関数に返り値を渡してdata.frame形式で出力してやれば見れる
365:１３２人目の素数さん
18/04/04 12:43:50.08 LbKgW3kd.net
>>353
>>355が言う直接的な参照
> fisher.test(matrix(c(1,120,130,2),2))$p.value
[1] 1.691912e-69
366:１３２人目の素数さん
18/04/04 13:06:17.36 PZp+DZN4.net
352です。よく分かりました！
ありがとうございます！
367:１３２人目の素数さん
18/04/04 14:44:33.05 LbKgW3kd.net
>>357
技術的な助言をしたけど、学術的に言えば、
p < 0.01 は全て p < 0.01 として、具体的なp値を考える意味はないと思うよ。
一部の例外的な研究分野を除いて(e.g., 遺伝統計学)。
368:１３２人目の素数さん
18/04/04 15:18:50.72 PZp+DZN4.net
はい、まさにその例外的な分野で使おうとしてます。ありがとうございます。
369:￥
18/04/07 06:52:38.98 yx+HETs3.net
￥
370:￥
18/04/07 06:53:00.80 yx+HETs3.net
￥
371:￥
18/04/07 06:53:20.77 yx+HETs3.net
￥
372:￥
18/04/07 06:53:40.58 yx+HETs3.net
￥
373:￥
18/04/07 06:54:01.92 yx+HETs3.net
￥
374:￥
18/04/07 06:54:24.32 yx+HETs3.net
￥
375:￥
18/04/07 06:54:47.76 yx+HETs3.net
￥
376:￥
18/04/07 06:55:08.41 yx+HETs3.net
￥
377:￥
18/04/07 06:55:30.46 yx+HETs3.net
￥
378:￥
18/04/07 06:55:54.53 yx+HETs3.net
￥
379:１３２人目の素数さん
18/04/07 10:24:54.55 maRES4NF.net
>>358
多重比較だと意味あるかも
380:１３２人目の素数さん
18/04/19 21:35:35.54 GVMUXyX9.net
Rのガンマ関数はいくつでオーバーフローするかやってみた。
> i=1
> while(gamma(i)!=Inf){
+ i=i+1
+ }
Warning message:
In gamma(i) : value out of range in 'gammafn'
> i
[1] 172
> gamma(172)
[1] Inf
Warning message:
value out of range in 'gammafn'
> gamma(171)
[1] 7.257416e+306
381:１３２人目の素数さん
18/04/29 00:21:51.67 5dW+xNwa.net
matplot()で折れ線グラフ描いたときに、X軸をカテゴリで示したいのですが、
可能でしょうか？
例えばtemp <- c("0時間","8時間","24時間","48時間")として、
matplot()の引数にtempをとるやり方です。
他にもやり方あれば教えてください。
382:１３２人目の素数さん
18/04/29 01:55:59.75 PngyHOQZ.net
>>372
matplot(..., xaxt="n")
axis(1, at=seq(along=temp), lab=temp)
383:１３２人目の素数さん
18/04/30 16:50:51.65 t3vhzyao.net
>>373
遅くなりましたがありがとうございました。
できました！
384:￥
18/04/30 23:55:02.74 y1TqbRSE.net
￥
385:￥
18/04/30 23:55:22.50 y1TqbRSE.net
￥
386:￥
18/04/30 23:55:43.19 y1TqbRSE.net
￥
387:￥
18/04/30 23:56:02.89 y1TqbRSE.net
￥
388:￥
18/04/30 23:56:23.17 y1TqbRSE.net
￥
389:￥
18/04/30 23:56:40.16 y1TqbRSE.net
￥
390:￥
18/04/30 23:56:53.49 y1TqbRSE.net
￥
391:￥
18/04/30 23:57:13.43 y1TqbRSE.net
￥
392:￥
18/04/30 23:57:32.05 y1TqbRSE.net
￥
393:￥
18/04/30 23:57:53.02 y1TqbRSE.net
￥
394:１３２人目の素数さん
18/05/01 18:57:32.34 iUBwAKWd.net
特定の長方形の中に複数の長方形を最小面積で敷き詰める平面充填に関するパッケージってありませんかね
395:１３２人目の素数さん
18/05/01 19:35:44.26 iUBwAKWd.net
やっぱ自分でアルゴリズムを書いてみます
396:１３２人目の素数さん
18/05/02 00:09:39.04 xpXz1bJY.net
Package`deldir'
これどうなの？
397:１３２人目の素数さん
18/05/06 22:15:43.42 BK1CxH7U.net
# jonckheereテストを書いてみた
jonckheere <- function(L,
alternative = c("two.sided", "increasing", "decreasing"),
cat=TRUE){
# L : list of vectors A1,A2,...,Ak, with assumed tendency
How.Many.Greater.Pairs <- function(A,B){ # How many pairs of A[i] > B[j], count as 0.5 when equal,
n.a = length(A)
n.b = length(B)
how.many.greater.pairs = 0
for(i in 1:n.a){
for(j in 1:n.b){
how.many.greater.pairs = how.many.greater.pairs+ifelse(A[i]==B[j],0.5,A[i]>B[j])
}
}
return(how.many.greater.pairs)
}
Sum.of.Greater.Pairs <- function(L){ #L=list(A1,,,,Ak),A1 < A2 < A3,..,< Ak : vector
k = length(L)
comb = combn(1:k,2) # possible combinaition of pairs to compare
n.comb = ncol(comb) # how many combinations
J = 0 # sum of greater pairs
for(i in 1:n.comb){
J = J + How.Many.Greater.Pairs(L[[comb[1,i]]],L[[comb[2,i]]])
}
return(J)
}
J = Sum.of.Greater.Pairs(L)
n = sapply(L,length)
N = sum(n)
EJ = (N^2-sum(n^2))/4
VJ = (N^2*(2*N+3)-sum(n^2*(2*n+3)))/72
Z = (J-EJ)/sqrt(VJ)
alternative = match.arg(alternative)
p.value = switch(alternative, 'two.sided' = 2 * min(pnorm(Z), pnorm(-Z), 0.5),
'increasing' = pnorm(Z),
'decreasing' = pnorm(-Z))
if(cat){
cat( 'p.value = ', p.value,'\n')
cat('alternative hypothesis: ' ,alternative,'\n')
}
invisible(p.value)
}
398:１３２人目の素数さん
18/05/22 20:34:14.53 8iM7y57s.net
代入は「<-」派？「=」派？
399:１３２人目の素数さん
18/05/22 21:16:34.39 iB1pjrmI.net
>>389
実態調査か何か？
<-と=は挙動が違う場合があるので、使い分けていますが、
代入はどっちか�
400:ﾆ問われたら、無論 <- または -> なお、 > 1 -> x これはエラーにならないけど、 > 1 = x 1 = x でエラー: 代入の左辺が不正 (do_set) です これはエラー
401:１３２人目の素数さん
18/05/23 11:17:56.08 OSJ/4EBd.net
>>390
俺は基本=派。
関数の定義は
z.test <- function(x,n=16,sigma=1){
z=sqrt(n)*mean(x)/sigma
2*pnorm(abs(z),lower=FALSE)
}
と書いている。
402:１３２人目の素数さん
18/05/23 11:24:57.67 MGQGuwX9.net
＞1 = x でエラー
当たり前
そんな使い方なんてするかよ
他言語と同じく=一文字の方がすっきりしてイイ
403:１３２人目の素数さん
18/05/24 14:28:38.60 ExPgBsbL.net
こういうのが紛らわしいから、俺は = 推奨。
x <- 1
if(x <- 1) print('YES')
if(x < -1) print('YES')
404:１３２人目の素数さん
18/05/24 14:36:00.11 EQ5K0CF7.net
だよねぇ
<-良くない
405:１３２人目の素数さん
18/05/25 08:05:58.04 ZHt2t+40.net
やったことなかったので関数の初期値設定に<-を使うとどうなるかやってみた。
まず、= の場合
> 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
> }
406:１３２人目の素数さん
18/05/25 08:11:48.02 ZHt2t+40.net
俺は見栄えがいいと思って関数定義には<-を使っているけど
= でも通常に動作する。
> 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
<- 推奨の人に聞きたいのだけど
<- でないと動作しないってことあるのだろうか？
407:１３２人目の素数さん
18/05/25 09:41:36.92 GCdXSt8a.net
<-は、打つのが
=より面倒
うっとおしい
408:１３２人目の素数さん
18/05/25 20:19:41.68 dXqzteX1.net
RStudioつかってりゃ［Alt+-］で簡単入力
それと引数の指定は代入じゃないと思うんだが感覚が違うよかな？
409:１３２人目の素数さん
18/05/25 21:33:31.44 QNMt6z2O.net
>>397
essだとアンダースコアを入れるを勝手に「<-」になる。
>>395
それは代入ではなく、関数の規定値の設定。
規定値の設定は「=」と決められているので、「<-」はアウト。
ただし、関数を実行するときには「<-」を使うことができる。
> mean(x<-1:10)
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10
410:１３２人目の素数さん
18/05/25 22:26:21.26 Cz/+g7h9.net
=に拘り続けて、<-に変える瞬間の君たちがたまらないよ
411:１３２人目の素数さん
18/05/28 16:26:17.30 d/09kgU6.net
>>399
ありがとうございます。
そんな使い方ができたのですね。
こんなのしか知りませんでした。
> f <- function() {
+ x<<-1:10
+ mean (x)
+ }
> f()
[1] 5.5
> x
[1] 1 2 3 4 5 6 7 8 9 10
412:１３２人目の素数さん
18/05/28 20:56:51.59 Osxttqv4.net
一画面しかグラフがないのに
Hit <Return> to see next plot:
と出るの鬱陶しいな、と思っていたら
par(ask=FALSE)
と設定しておけばいいんだな。
413:１３２人目の素数さん
18/05/29 18:08:52.29 ZCWRCH5Y.net
同名のdataがあるのでパッケージを::で指定してもうまくいかなかった。
data(netmeta::parkinson)
でなくて
data(parkinson, package = 'netmeta')
とするんだな。
414:１３２人目の素数さん
18/05/30 06:16:03.29 fLd3NENr.net
hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
だと、1色だけど
hist(rnorm(100),col=sample(colours(),(sample(1:10,1))))
なら1～10色で表示される。
415:１３２人目の素数さん
2018/05/31
416:(木) 21:37:03.79 ID:BfVjjX7C.net
417:１３２人目の素数さん
18/05/31 23:16:15.39 grs1zCKo.net
> hist(rnorm(100),col=rgb(runif(1),runif(1),runif(1),runif(1)))
何色に増やしてもいいがな。
何をやりたいのかな？
hist(rnorm(5000),col=apply(matrix(runif(80),4), 2, function(x){rgb(x[1],x[2],x[3],x[4])}))
418:１３２人目の素数さん
18/06/01 21:54:20.67 Ef899k/0.net
:::てどういう時に使うんだろ？
ソース読みたくて
library(BayesFactor)
ttestBF_indepSample
では表示されなかったが、
library(BayesFactor)
BayesFactor:::ttestBF_indepSample
だと出てきた。
419:１３２人目の素数さん
18/06/04 13:04:13.10 3QOAvZa7.net
Stanも覚えなきゃいけないのかよーくそがー
420:１３２人目の素数さん
18/06/04 13:05:22.57 3QOAvZa7.net
::なら名前空間カブールよけよけの術だよね
421:１３２人目の素数さん
18/06/04 14:24:00.43 d7l4LqX9.net
>>408
とりあえずこれ
URLﾘﾝｸ(github.com)
422:１３２人目の素数さん
18/06/05 12:44:14.88 LfwL/Kok.net
>>409
:::と3個なんだよなぁ
423:１３２人目の素数さん
18/06/05 13:06:03.27 CHJNgrVq.net
>>411
グダグダ言っていないで、さっさとヘルプを参照すればいいじゃないか？
?':::'
424:１３２人目の素数さん
18/06/05 19:45:28.11 LfwL/Kok.net
>>412
ありがとう。
stats:::t.test.default
でｔ検定のソースが見られた
425:１３２人目の素数さん
18/06/11 00:14:18.19 K1zGYlRd.net
関数を値で上書きしてしまった場合に元に戻す方法を教えて頂きたく。
例)
rm <- 1
としてしまった場合にrm関数を元に戻したいです。
よろしくお願いします。
426:１３２人目の素数さん
18/06/11 07:21:35.69 +ayC4RGv.net
>>414
直後ならRstudioでctrl+zで元に戻す。
再起動してたらこの方法は無理だろな
427:１３２人目の素数さん
18/06/11 08:24:34.24 2ZHCUqeW.net
>>414
自分で作ったのではない関数の rm だったら rm(rm) で数値ベクトルオブジェクトが消えて、関数の rm が見えるようになる。
428:１３２人目の素数さん
18/06/11 15:01:37.34 mlFyU0v4.net
>>414
> 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()で大丈夫だろう
429:１３２人目の素数さん
18/06/11 17:04:35.54 Z+okZT62.net
>>414
既に説明されてるけど上書きじゃなくて別オブジェクト扱いになるからrmで消すか明示的にパッケージを::演算子で関数を指定してやればおけ
自作関数だと上書きされちゃうけど
430:１３２人目の素数さん
18/06/11 19:33:14.14 rfEFkwvW.net
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)
431:１３２人目の素数さん
18/06/11 19:51:41.60 Z+okZT62.net
同じ名前空間にあるなら上書きされちゃうよ
.Gloalenv（だったかな？）に自作関数があり、そこに上書きして変数にしてるから消したら消えるだけで戻せない
432:１３２人目の素数さん
18/06/11 19:55:10.50 Z+okZT62.net
ちょっと古いけど、この辺りを読むとなぜそうなるかが分かると思う
（Rの）環境問題について その１。
URLﾘﾝｸ(qiita.com)
433:１３２人目の素数さん
18/06/11 20:04:41.37 mlFyU0v4.net
>>419
戻せます。
> 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),
[以下略]
434:１３２人目の素数さん
18/06/11 20:06:37.76 mlFyU0v4.net
>>420
>.Gloalenv
case sensitiveなRでそれはないのでは？
> help(".GlobalEnv")
435:１３２人目の素数さん
18/06/11 20:19:14.99 Z+okZT62.net
フォローありがとう
436:413
18/06/11 21:00:49.78 K1zGYlRd.net
>>415-423
みなさまご返信ありがとうございました。
例示したものがbase::rmだったので自作関数のことまでは想像だにしませんでした…(がとても良い勉強になりました。)
自環境でも上書きした関数を元に戻すことが出来ることを確認しました。
437:１３２人目の素数さん
18/06/12 07:42:40.88 9RfEtlLW.net
>>422
戻したいのが自作関数
rm = function(x) asin(sqrt(x))
だったら上書きされて無理じゃないの？
ゴミ箱から消去したファイルを復活させるソフトもあるから
PCに詳しければ可能かもしれないけれど
Rの機能としては復活させるのは無理と思う。
438:１３２人目の素数さん
18/06/12 11:00:45.94 cD/mQB+6.net
自作のrm関数はbaseパッケージとは別に定義されるから
base::rm(rm)で自作の方のrmは消えてbase::rmが「復活」します。
439:１３２人目の素数さん
18/06/12 12:07:52.53 cD/mQB+6.net
ああすまん
自作の方を上書きした場合は、そっちの復活は無理ですね。
440:１３２人目の素数さん
18/06/12 14:36:07.63 VF7OcVBc.net
>>426
ちょっと意味がわからない。
自作の関数なら定義を再実行するだけなのに、
消える消えないって何が問題なんだ？
441:１３２人目の素数さん
18/06/13 14:39:27.58 ZXkz+qgi.net
>>425
関数だけでなく既定値でも同じ。
> pi
[1] 3.141593
> pi=2
> pi
[1] 2
> rm(pi)
> pi
[1] 3.14159
442:１３２人目の素数さん
18/06/14 09:27:55.98 mxBGyFKT.net
共同ツール 1
URLﾘﾝｸ(seleck.cc)
URLﾘﾝｸ(trello.com)
ボードのメニュー → Power-Upsから拡張可能 Slack DropBoxなど
Trello Chrome拡張機能 elegant
URLﾘﾝｸ(www.kikakulabo.com)
trelloのオープンソースあり
共同ツール 2
URLﾘﾝｸ(www.google.com)
共同ツール 3
URLﾘﾝｸ(slack.com)
URLﾘﾝｸ(www.dropbox.com)
URLﾘﾝｸ(bitbucket.org)
URLﾘﾝｸ(ja.atlassian.com)
URLﾘﾝｸ(www.sketchapp.com)
URLﾘﾝｸ(photoshopvip.net)
URLﾘﾝｸ(goodpatch.com)
Trello Chrome拡張機能プラグイン集
URLﾘﾝｸ(chrome.google.com)
Slackプラグイン集
URLﾘﾝｸ(slack.com)
Sketchプラグイン集
URLﾘﾝｸ(sketchapp.com)
URLﾘﾝｸ(supernova.studio)
443:１３２人目の素数さん
18/06/16 19:29:24.89 10mr4q2e.net
はっとデシャンｸﾞ
444:１３２人目の素数さん
18/06/17 19:06:33.58 OYjqtCQI.net
>>307
知恵袋に漸化式があったので計算スクリプトを書いてみた。
表がでる確率ｐのコインを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
445:(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] }
446:１３２人目の素数さん
18/06/18 09:16:00.35 JdABdagZ.net
>>307
この例で、答えとして何が返ってきて欲しいのか？
2で合ってる？
447:１３２人目の素数さん
18/06/18 21:23:58.67 fHaLluDH.net
>>434
1が5回以上出現しているからTRUE（数字なら１）が返り値。
448:１３２人目の素数さん
18/06/19 10:20:04.23 fcA67BpN.net
>>435
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があるか?
}
449:１３２人目の素数さん
18/06/19 21:57:17.85 QuKrVrE8.net
>>436
whichとdiffを用いての解法のお返事ありがとうございます。
実行時間を>303のスクリプトと比べてみました。
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の使い方が勉強になりました。
時間を割いてスクリプトを作成していただいてありがとうございます。<(_ _)>
450:１３２人目の素数さん
18/06/21 10:12:41.96 Ze2kEGUX.net
Rstudioで日本語の入力ができないのですが、解決する方法知ってる方いませんか?
具体的には、IME自体をオンにすることができません。
日本語の表示やコピペは問題なくできます。
環境
Linux - Utubntu
RStudio Version 1.0.143
451:１３２人目の素数さん
18/06/21 10:26:53.84 Fds/4hNR.net
RStudioのQtが古いのでパッチ当てないと日本語入力できないのが現状ですわ
パッチはUbuntu 16.04.LTS用にしかないので他のバージョンならRStudioサーバー使った方が早いかと
なお、パッチはRStudio 日本語とかでグクってみて
452:１３２人目の素数さん
18/06/23 17:38:43.39 aZCZP6wm.net
windows版のRStudioでもたまに日本語入力ができなくなる。
別のソフトでIMEのON/OFFをやってから戻ると直るけど。
453:１３２人目の素数さん
18/06/23 20:18:55.21 4XxgOFIA.net
>>440
うちも同じ。
454:１３２人目の素数さん
18/06/23 21:30:47.03 JcH4CpHc.net
Windowsはストアアプリでも同じ症状でるから持病なんじゃないかと思う
Ubuntuはfcitxのパッチ当てとけばWindowsのような症状は出ないので快適
455:１３２人目の素数さん
18/06/25 13:23:09.19 rmZ6BvQE.net
これをやってみた。
URLﾘﾝｸ(img-cdn.jg.jugem.jp)
スクリプトは
URLﾘﾝｸ(imagizer.imageshack.com)
456:0/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
457:１３２人目の素数さん
18/06/26 07:08:45.35 Na/Ih9Bj.net
問題
９９人の囚人がいます。彼らの頭に１～１００までのナンバーカードが貼りつけられた帽子をランダムにかぶせます。
他人の帽子は見ることができても、自分の帽子は見ることができません。
帽子の数は全部で１００なので、一つ使われずに余ります。
そのナンバーは囚人達にはわからないようにしておきます。
この状況で、囚人たちに一斉に自分のナンバーを宣言させて、全員が正解だったら釈放するという賭けをします。
囚人たちには帽子をかぶせられる前に相談タイムが設けられています。
どういう戦略を取れば、助かる確率を最も高くできるでしょうか？
URLﾘﾝｸ(000013.blogspot.com)
解答を読んでも数理が理解できないが
確率が0.5になるのはシミュレーションできた。
# URLﾘﾝｸ(000013.blogspot.com)
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()))
458:１３２人目の素数さん
18/07/11 16:20:29.88 zFHa28EV.net
>>438
遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
詳しくはUbuntu Weekly Recipeの第527回を読んでくだされ
459:１３２人目の素数さん
18/07/22 11:22:38.35 Ott8rTSz.net
覆面算 RED + WHITE = COLOR
どうもＣでやったらうまくいかないのでＲでやってみた。
# 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
460: 0 2 0 5
461:１３２人目の素数さん
18/07/22 12:50:02.38 Ott8rTSz.net
# AAB
# × CD
# ------
# CCE
# BFAA
# ------
# BEBDE
f <- function(A,B,C,D,E,F){
(A*100+A*10+B)*D==C*100+C*10+E &(A*100+A*10+B)*C==B*1000+F*100+A*10+A &
(C*100+C*10+E)+(B*1000+F*100+A*10+A)*10==B*10000+E*1000+B*100+D*10+E
}
library(gtools)
perm=permutations(n=10,r=6,v=0:9)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=f(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6])
}
indx=which(re==TRUE)
perm[indx,]
462:１３２人目の素数さん
18/07/23 11:09:01.32 SDW/E7TY.net
大変初歩的な質問で恐縮ですが、ご教示いただければ幸いです。
PERT分布を学んでいるのですが、ベータ分布のスケーリングの仕方が分かりません。
最小値a、最頻値ｂ、最大値ｃ として、a, b, c から求めたパラメータをa1、a2 とします。
このとき、
dbeta(x, a1, a2)*(c-a)+a
とすると、当然ながら横軸の値の範囲は0から1で、縦軸の値のみ大きくなってしまいます。
横軸の値の範囲がaからcで、縦軸の値は確率密度のままにするには、どうしたらよいのでしょうか？
463:１３２人目の素数さん
18/07/23 12:15:25.71 llIMc79p.net
え？aだけ右に動かしたいならこれでいいのでは？
dbeta(x-a, a1, a2)
464:１３２人目の素数さん
18/07/23 18:17:20.65 SDW/E7TY.net
>>449
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の一様分布のようなグラフになってしまいます。
465:１３２人目の素数さん
18/07/23 20:05:23.59 5MGL8f5L.net
4パラべーたへの変換は下記かな？
URLﾘﾝｸ(en.wikipedia.org)
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))
466:１３２人目の素数さん
18/07/23 21:05:56.52 +HoGXx9d.net
>>451
447です。
450様、まさにこれ、これです。どうもありがとうございました。
おかげさまで大変助かりました。
448様にも、ご助言下さいましたこｔ、厚く御礼申し上げます。
467:１３２人目の素数さん
18/07/27 15:50:12.19 E7dz8jqT.net
lmtestのパッケージをインストールしようとすると不正なマルチバイト文字がありますと出てくるのですがどうすればいいですか？
468:１３２人目の素数さん
18/07/27 17:33:01.36 ufA3ZDTu.net
>>453
Windows 10 （Rの日本語ヘルプはインストールしていない）+ R studioだけど、そのメッセージ出ないな。
469:１３２人目の素数さん
18/07/27 23:29:04.58 klPq25x/.net
ユーザー名が日本語かな
470:１３２人目の素数さん
18/07/28 12:56:05.38 lhV9Awbf.net
職場のSPSS信者たちにRを認めてもらうにはどうしたら良いだろう
「安かろう悪かろう」「俺が知らないものは三流」的な考え方の人が多くて肩身が狭い…
471:１３２人目の素数さん
18/07/28 17:53:14.05 XHIBT/Gx.net
>>456
>>453のようなエラーが出たら「無料のものはやっぱりダメね�
472:v 有償(SPSS) →信頼できる ボランティア→いい加減 と言う思考らしい
473:１３２人目の素数さん
18/07/28 18:33:45.13 l8Bhwsl9.net
>>455
ユーザー名はアルファベット1文字です
studioインストールしてないんですけど、しないとダメですか？
474:１３２人目の素数さん
18/07/28 21:48:13.01 hNHqu0EH.net
他のパッケージはインストールできるの？
475:１３２人目の素数さん
18/07/29 07:46:06.57 u2P99O/7.net
linuxなら環境変数LANGをCにすれば解決する問題だからwindowsでも似たような対処すれば行けそう。
476:１３２人目の素数さん
18/07/29 19:38:42.16 hjP1s1BE.net
これかな？
> Sys.setlocale(locale="C") # locale の変更（USの設定にしないと表示が変になる）
[1] "C"
元に戻すには
> Sys.setlocale(locale="Japanese_Japan.932")
477:１３２人目の素数さん
18/08/02 20:59:18.66 DvtHDake.net
すみません、教えて下さい。
フィッシャーの正確確率検定で、
適合度（一様性）の検定ってできますか？
できるのは独立性の検定のみ？
Rでいろいろ試したみたんですが、
2群ないとエラーになる。。
478:１３２人目の素数さん
18/08/02 21:16:50.41 GWs8Hfcd.net
適合度ならカイ二乗検定でやればいいのでは？
479:１３２人目の素数さん
18/08/02 21:52:11.65 DvtHDake.net
>>463
レス感謝です。
そうなんですが、
期待度数が5未満になってしまった場合は、
カイ二乗検定は向かないんですよね？
480:１３２人目の素数さん
18/08/02 21:54:12.25 82PNyNf0.net
それで他の方法を探しています
481:１３２人目の素数さん
18/08/02 23:02:57.38 clch7kxN.net
>>464
こういうので
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
482:１３２人目の素数さん
18/08/03 09:04:31.58 VGe+pklE.net
>>464
適合度を見るのに関係ないだろ。
> 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, simulat
483:e.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未満の数があっても関係がない
484:１３２人目の素数さん
18/08/03 10:29:47.68 jRZN4rrt.net
>>467
警告が出るのは独立性検定の場合でしたね。
485:１３２人目の素数さん
18/08/04 02:59:15.87 8dNre62F.net
463です。
カイ二乗検定で、期待度が5未満が多いのは不適切なのかと思い込んでおりましたが、
それは独立性の検定の場合のことであり、
適合度（一様性）の検定の場合は、
気にしなくて良かったのか。。
ならカイ二乗検定でやります。
いろいろ教えて下さり、ありがとうございます。
486:１３２人目の素数さん
18/08/04 19:49:50.94 3IVOmTrE.net
なんでお前らSPSS使わないの？
URLﾘﾝｸ(www.cdleidyba.lt)
487:１３２人目の素数さん
18/08/04 20:37:20.39 kYh4t8ZN.net
>>470
Stanとかに対応してんの？
488:１３２人目の素数さん
18/08/04 22:03:56.61 kTWj7ufe.net
>>470
化石みたいなバージョンやな
489:１３２人目の素数さん
18/08/05 00:56:32.45 biDT5wEJ.net
>>470
いまSPSSを使う理由って何なの？煽りじゃなくマジで
490:１３２人目の素数さん
18/08/05 14:27:54.32 dSJLM/da.net
maple使えよ
491:１３２人目の素数さん
18/08/05 19:10:34.43 mULz5b3r.net
>>473
周囲のspssユーザを見ると「指導教員がspssしか使えないので、
自分もspssしか使えない、
新しい（学術系）ソフトは誰かに丁寧に根気よく指導してみらえないと無理」、
これの再生産。
492:１３２人目の素数さん
18/08/05 19:11:30.48 zaJ/WJYb.net
>>470
研究室に予算がないから
493:１３２人目の素数さん
18/08/05 19:12:42.57 mULz5b3r.net
>>476
予算がなかったら、Rを使え
494:１３２人目の素数さん
18/08/05 19:13:59.36 mULz5b3r.net
>>470と>>476を間違えたorz
495:１３２人目の素数さん
18/08/05 19:33:14.43 mTZcaJ/n.net
予算がないならケーキを食べればいいのに
496:１３２人目の素数さん
18/08/08 13:11:20.17 7/3l/PxD.net
再帰呼び出しに関数名の代わりにRecall使うと処理に要する時間が増えることに気づいた。
フィボナッチ数列の再帰呼び出し
# 関数名での呼び出し
fibo <- function(n){
if(n==1|n==2) return(1)
else fibo(n-1)+fibo(n-2)
}
# Recallを使っての呼び出し
fiboR <- function(n){
if(n==1|n==2) return(1)
else Recall(n-1)+Recall(n-2)
}
> system.time(fibo(30))
user system elapsed
4.27 0.02 4.36
> system.time(fiboR(30))
user system elapsed
6.65 0.00 6.70
497:１３２人目の素数さん
18/08/09 17:47:22.65 l1M8GWWv.net
cumsumを再帰関数で書いてみた。
何度も試行錯誤
# cumsum with recursive call
cumsumR <- function(v){
res=numeric(length(v))
cumsumR_sub <- function(v,res,i){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
cumsumR_sub(v,res,1)
}
cumsumR(1:10)
とりえあず、動作した。
> cumsumR(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
498:１３２人目の素数さん
18/08/09 18:31:24.68 l1M8GWWv.net
簡略化できた
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
>
499:１３２人目の素数さん
18/08/10 18:41:42.48 Hlm8Oe3x.net
ifelseの動作がどうも納得できないでご教示いただきたいのですが。
これってifelseの仕様でしょうか？バグでしょうか？
> x=2:1
> if(x[1]<=x[2]) x else x[2:1]
[1] 1 2
> ifelse(x[1]<=x[2],x,x[2:1])
[1] 1
500:１３２人目の素数さん
18/08/10 18:55:39.89 vdkgfABT.net
>>483
ifelse() はベクトル演算できる関数。
だから、戻り値の要素数は、test部分の要素数に一致する。
testが1つなので、x[2:1]の1つ目の要素が返される。
501:483
18/08/10 18:58:51.59 vdkgfABT.net
testの要素を2つにすれば、出力も2つになる。
> 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
502:１３２人目の素数さん
18/08/10 23:06:17.77 Hlm8Oe3x.net
>>484
ありがとうございました。
バグではなくてそういう仕様だったのですね。
こういう使い方ができるということがわかって勉強になりました。
> x=2:1
> ifelse(c(x[1]!=x[2],x[1]==x[2]),1:2,3:4)
[1] 1 4
503:１３２人目の素数さん
18/08/10 23:17:14.10 VlOUWHrC.net
教わったので早速、
ifelse() はベクトル演算できると、再利用されるの動作確認。
> ifelse(c(TRUE,TRUE,FALSE,FALSE,TRUE,FALSE,FALSE),1:2,11:15)
[1] 1 2 13 14 1 11 12
504:１３２人目の素数さん
18/08/10 23:45:24.37 VlOUWHrC.net
２種類の文字の変換もifelseでできるんだなぁ。
> mat
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 0
[3,] 1 1 0
[4,] 1 1 0
[5,] 1 1 0
[6,] 1 1 0
> print(apply(mat,2, function(x) ifelse(x==1,'真','偽')),quote=F)
[,1] [,2] [,3]
[1,] 真 真 真
[2,] 真 真 偽
[3,] 真 真 偽
[4,] 真 真 偽
[5,] 真 真 偽
[6,] 真 真 偽
505:１３２人目の素数さん
18/08/11 00:28:09.90 OesSEWnz.net
行列も次元付きのベクトルだからもっと簡略化できた。
> (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
>
506:１３２人目の素数さん
18/08/11 18:04:43.46 05wxJPem.net
RやSPSSって名前の試験まであるのね
507:１３２人目の素数さん
18/08/16 23:26:21.81 3mKQ9SP5+
>>469
適合度検定の場合は、不適の基準は期待度数5未満ではなく、2.5未満では?
つまり一様性の検定なら、全ての期待度数が2.5以上なら良いのでは?