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以上なら良いのでは?
508:１３２人目の素数さん
18/08/22 14:02:15.21 rznk0lAS.net
規約分数にするパッケージが探せなかったので自分でスクリプトを組んでみた。
エラー処置は省略。
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
509:１３２人目の素数さん
18/08/23 01:12:39.77 0Wz4IoKN.net
rvest?でログイン必要なとこをスクレイピングする時って
login_page <- html_session("URLﾘﾝｸ(xxxxx"))
login_form <- html_form(login_page)[[1]] %>%
set_values(AAA="xxxx", BBB="xxxx")
というので行けるはずなんですが
最終行のAAAとBBBってそれぞれhtmlタグのname=""のとこからとってくるんですよね?
これにここがname="user[email]"みたいに[]という記号はいってたらどうすればいいでしょうか?
510:１３２人目の素数さん
18/08/25 17:48:59.21 MdsDkupV.net
>>493
よろしくおねがいいしまーーーーーーーす
511:１３２人目の素数さん
18/08/26 22:01:57.06 2bj/HrEX.net
>>494
質問の意味がいまいち分からないから、誰も助言できないのでは？
> x <- 'user[email]'
> x
[1] "user[email]"
rvestパッケージを使ったことがないけど、
記号が入っていたらなぜ問題があるのかよく分からない。
512:１３２人目の素数さん
18/08/26 22:53:43.10 WnqFRbi9.net
というかどこのページのスクレイピングしたいか晒して
仮でもいいから
513:１３２人目の素数さん
18/08/27 06:16:23.20 x8E4FG0O.net
>>495 ちがいます AAAがどれかよくみてみてください ””はありません
515:１３２人目の素数さん
18/08/27 06:19:43.00 x8E4FG0O.net
ちょっと直接はることはできないのでタグだけはります。
<input class="form-control" autofocus="autofocus" placeholder="Email address" required="required" type="text" value="" name="user[email]" id="user_email">>>496
516:１３２人目の素数さん
18/08/27 13:25:39.67 SdOxff6m.net
>>497
もっと分かるように説明しないと、追試できる情報も提供しないし。
もしかして変数名などに[]が入っている場合 にどうしたらよいかってこと？
> `user[email]` <- 1:5
> `user[email]`
[1] 1 2 3 4 5
517:１３２人目の素数さん
18/08/27 15:43:37.46 x8E4FG0O.net
>>499
多分行けました
超感謝!!
518:１３２人目の素数さん
18/08/27 20:20:32.13 5a+trmgU.net
ある問題のシミュレーションしようと思って
問題文の記号のまま
q <- function(x)
....
とやって
q(100)と入力すると
Rが終了することに気づいた。
519:１３２人目の素数さん
18/08/27 22:56:54.68 a+8R0Tu1.net
base::q()
を実行すると爆弾出るね
520:１３２人目の素数さん
18/08/29 22:45:10.53 BcwFyR33.net
ちょっとした疑問です。
空ベクトルの検出って長さ０以外に検出方法ってあるでしょうか？
> 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
521:１３２人目の素数さん
18/08/31 22:46:15.16 IWQvY6FL.net
４点の座標を入力するとそれらを結ぶ四面体の体積を求めるスクリプトを書いてみた。
高さはパッケージ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
多分、正常に動作していると思う。
522:１３２人目の素数さん
18/09/01 00:25:33.05 52Ub52jp.net
>>504
行列式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
523:１３２人目の素数さん
18/09/01 01:41:47.56 qG52f2Ee.net
ベクトルの三重積を教わったので、パッケージ pracma の外積crossを使った
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)
}
４行で済んだ。
>>505
ありがとうございました。
パッケージに頼らずに計算できたのですね。
524:１３２人目の素数さん
18/09/03 18:55:17.77 S47YTHgP.net
データを解析する前にさらっと特徴を見たい時、皆さんはどんなコマンドを使っていますか？
私が思いつくのは
summary
boxplot
hist
pairs
です。こんなのも良いよってのがあったら教えてくださいm(_ _)m
※ライブラリの使用有無は問いません
525:１３２人目の素数さん
18/09/03 19:38:36.45 fowZfPON.net
>>507
BESTパッケージのplotPost
histに加えて95%CIを表示してくれる。
526:１３２人目の素数さん
18/09/03 20:05:14.66 OGj8hrn2.net
>>507
出力をデータフレーム型にしたいときはskimrパッケージ
527:１３２人目の素数さん
18/09/03 20:07:16.33 OGj8hrn2.net
あと、まだ試してないけsummarytoolsパッケージも面白そう
528:１３２人目の素数さん
18/09/09 22:35:22.23 9XY+z1xx.net
>>503
良い方法が見つからなかったので !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
529:１３２人目の素数さん
18/09/09 23:21:54.06 1udV8jVZ.net
>>511
なぜ if (length(x) == 0) と書かないのか？
可読性って知ってるか？
530:１３２人目の素数さん
18/09/10 16:54:23.64 89eIPezd.net
>>512
文脈読めば
!length(x)で空白ベクトル
の流れ
531:１３２人目の素数さん
18/09/10 18:22:05.74 gLTPxqtw.net
>>513
length は数値を返す関数なのに、なぜわざわざ論理演算をするんだ？
プログラムは、ここの議論を知らない人が読む可能性のほうが高いんだぞ。
532:１３２人目の素数さん
18/09/10 19:12:52.95 ZYY4OYkH.net
>>514
正規分布が一様分布より大きい期待値の算出に
mean(rnorm(100)>runif(100))
って書かない？
俺は
sum(ifelse(rnorm(100)>runif(100),1,0))/100
って書いたりしないけど。
533:１３２人目の素数さん
18/09/10 19:43:25.60 SKQ8XoAj.net
>>514
え、常套手段やん？！
534:１３２人目の素数さん
18/09/10 19:50:35.36 ZYY4OYkH.net
>>516
論理値を数値に置き換えて計算しているから
数値を論理値にしても別に違和感がない。
可読性は慣れの問題。
535:１３２人目の素数さん
18/09/10 19:54:17.32 SKQ8XoAj.net
>>517
アンカー間違えとる？
だから、そんなの常套手段やん？
536:１３２人目の素数さん
18/09/10 21:30:23.37 gLTPxqtw.net
>>515
それとこれとは話が別だ。
logical は TRUE か FALSE の二値しかとらないから、それから数値への変換は自明。
多種の値をとる数値に論理演算をするのはいただけない。
537:１３２人目の素数さん
18/09/10 21:36:20.92 gLTPxqtw.net
>>517
is.empty.vector のような関数があるならそれでよいが、length を使っているのだから、それを勝手にベクトルが空かどうかの判断に使っているのはあなたであって、他人はそのようには考えない、といっているのだ。
自分だけしかこーどを見ないならべつに構わないが、こういうところに晒すのはよくない。
ましてや常套手段などと言って他人に教え込むのはやめてもらいたい。
538:１３２人目の素数さん
18/09/10 21:51:19.40 5QS5/GHY.net
>>520
常套手段というのはmeanの話じゃね？
539:１３２人目の素数さん
18/09/10 21:55:50.32 5QS5/GHY.net
whileの中は１でも２でもよくね？
n=0
while(1){
if(n>10) return(10)
n=n+1
}
n=0
while(2){
if(n>10) return(10)
n=n+1
}
540:１３２人目の素数さん
18/09/10 22:00:57.10 5QS5/GHY.net
>511は
配列を逆順に並べる再帰呼出しのコード。
Rには
541:revという関数がある。 そのコードみてみ！ lengthを真偽判定に使ってますがな。 > base:::rev.default function (x) if (length(x)) x[length(x):1L] else x
542:１３２人目の素数さん
18/09/10 22:18:05.27 DJR6mPp+.net
>>520
>523みた？
if(length(x)!=0)になってないよ。
543:１３２人目の素数さん
18/09/11 00:19:39.45 WPUC3B84.net
>>520
>こういうところ
誰も鵜呑みにしない便所の落書き。。。
544:１３２人目の素数さん
18/09/11 08:48:42.34 QUqp/jpE.net
!を引数が数値のときは0か否かを返す関数と理解すれば
if(!0) print(!1)の結果もサクッとわかる。
545:１３２人目の素数さん
18/09/11 09:08:32.20 QUqp/jpE.net
数値を非零かどうか論理値に変換して処理するかは
関数によるな。
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)
>
546:１３２人目の素数さん
18/09/11 09:12:56.04 QUqp/jpE.net
!ってベクトル対応しているな
> which(!(-1:1))
[1] 2
> !(-2:2)
[1] FALSE FALSE TRUE FALSE FALSE
547:１３２人目の素数さん
18/09/15 08:05:49.52 Vl7XZ52q.net
!をnotではなくis.falseとかis.zero変換すれば( ・∀・)ｲｲ!!
548:１３２人目の素数さん
18/09/15 16:42:00.50 h0gUCZ3r.net
>>529
is.zero だと？
だったら最初から == 0 てよいじゃないか。
549:１３２人目の素数さん
18/09/15 16:45:38.59 h0gUCZ3r.net
>>527
だいたい、if や while は関数じゃないし。
while (1) などと書くのはCに毒されているんじゃねーの？
Rなら while (TRUE) のほうが良いし、repeat というのもある。
550:１３２人目の素数さん
18/09/15 16:45:41.31 Vl7XZ52q.net
>>530
スクリプトを読むときに脳内変換する癖をつけておけば
可読性が悪いと思わずにすむ。
551:１３２人目の素数さん
18/09/15 16:46:13.42 Vl7XZ52q.net
>>520
>523みた？
if(length(x)!=0)になってないよ。
552:１３２人目の素数さん
18/09/15 17:19:35.45 h0gUCZ3r.net
>>533
Rのライブラリ？の書き方なんてあてにならないよ。
そもそも、rev も 509 の reverse も演算の回数を必要最低限にするなら、条件は length(x)) >= 2 または > 1 だ。
まあ if (length(x)) のほうが速かったのかも知れないが、常にそうなるとは限らないのだから、早すぎる最適化は諸悪の根源だ。
553:１３２人目の素数さん
18/09/15 19:30:34.55 Vl7XZ52q.net
可読性の次は速度かよ。
自分の流儀と違う { の使い方でも文句いいそう。
554:１３２人目の素数さん
18/09/15 19:44:27.43 VibLIqgl.net
0,1を1000万個作って
!
!=
>
==
での真偽判定を各々１０回する時間を出してみた。
> 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
555:１３２人目の素数さん
18/09/15 19:57:45.72 VibLIqgl.net
!x と x!=0 で１０回やったが、
> 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
結果は再現された。
556:１３２人目の素数さん
18/09/15 20:06:37.73 VibLIqgl.net
lengthが絡むと
> 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
御指摘の通り、!の負け
557:１３２人目の素数さん
18/10/04 17:02:48.04 247Ted9r.net
>>445
＞遅レスだが、Ubuntu 18.04LTSならデフォルトのIBus-mozcで日本語入力ができるそうだ。
＞詳しくはUbuntu We
558:ekly Recipeの第527回を読んでくだされ Ubuntu 18.04.1 LTSにアップグレードしたけどできない いままでスクレイピングにつかってrvestのログインコマンドもエラー出るようになったし最悪だ。。。。
559:１３２人目の素数さん
18/10/04 20:04:49.68 W6fWyA5T.net
アップグレードってことは文字入力がFcitxのままになってない？Fcitxならパッチ要だよ。
Voyager 18.04LTS（xubuntu 18.04LTS）ではパッチなしでRStudioに日本語入力できた。
560:１３２人目の素数さん
18/10/26 12:58:59.43 8LsI2NoU.net
昨日、PCを新しくして最新のR入れたんだけど、
作業ディレクトリの固定化できなくね？
前まではショートカットのプロパティで作業フォルダ指定してたんだけど、
新しいR、何度やっても作業ディレクトリがデフォルトのドキュメントに戻りやがる
どうすればいいか、誰がご教示を。
561:１３２人目の素数さん
18/10/26 13:04:50.23 I9GOHEF6.net
>>541
プロパティでコマンドラインオプションを消せばよい。
--cd-ほにゃらら
というのが付いてるでしょ。
562:１３２人目の素数さん
18/10/26 13:13:30.11 8LsI2NoU.net
>>542
あ、なるほど
ありがとうございます、多謝
563:１３２人目の素数さん
18/10/26 23:05:34.11 5yNEfvRx.net
xtsの取り扱い方について質問失礼します。
以下のコードで、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
564:１３２人目の素数さん
18/10/26 23:12:48.52 Xg1/ChR5.net
失礼、>>544の下のコードはこっちです
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
565:１３２人目の素数さん
18/10/26 23:38:29.60 oEWqNGH8.net
>>545
idx <- seq(10, length(time), by=24)
result <- data[idx]
じゃダメなのか？
566:１３２人目の素数さん
18/10/27 01:10:34.48 NQVW6zPX.net
>>546
ああ、これがいいですね！お恥ずかしい…
ありがとうございます
567:１３２人目の素数さん
18/11/01 22:44:27.08 xCdOvDq8.net
巨大数を扱えるというふれこみのRmpfrって正確じゃないな。
５０の階乗を計算させてみた
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:
URLﾘﾝｸ(www.wolframalpha.com)
568: 30414093201713378043612608166064768844377641568960512000000000000
569:１３２人目の素数さん
18/11/02 23:29:01.53 fMme23mF.net
こまけぇこたぁいいんだよ!!
570:１３２人目の素数さん
18/11/03 09:32:49.16 NiwdVATo.net
R3.5.xから、Windows環境だと日本語パスまわりでエラーがでてどうにもならない
enc2native()しても駄目だし
とりあえずcairo_png()だけでも…
571:１３２人目の素数さん
18/11/03 10:31:12.42 pqZePX+I.net
一度つかって日本語周りの処理でエラー出てきたんで未だに3.4.4浸かってる
572:１３２人目の素数さん
18/11/04 21:55:16.95 lTCeMsqQ.net
統計の質問はこのスレでいいんだろうか？
ある医院に１時間あたり平均５人の患者が来院し、その人数の分布はポアソン分布にしたがうとする。
1時間あたりの平均診療人数は６人で、一人あたりの診療時間は指数分布に従うとする。
診察までの平均の待ち時間は何時間か？
このシミュレーション解はこれであってる？
N=1e6
λ=5
μ=6
sum(rpois(N,λ)*rexp(N,μ))/N
> sum(rpois(N,λ)*rexp(N,μ))/N
[1] 0.8331036
60かけて、分にすると
[1] 49.86565
573:１３２人目の素数さん
18/11/08 12:35:54.73 wKTjJ6Fa.net
>>307
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
逆に４倍時間がかかるようになった
574:１３２人目の素数さん
18/11/08 14:05:00.34 1h1GfW+d.net
今さらだけどそれってrle使っちゃだめなの？
575:１３２人目の素数さん
18/11/08 15:13:38.19 UNhl34kg.net
coinの表裏を1と0で表して、その累積和を取ったベクトルを用意して
そのベクトルの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
576:１３２人目の素数さん
18/11/08 22:10:55.78 I6htNxzg.net
あっ、コード中の5はnに置き換えてね
577:１３２人目の素数さん
18/11/08 23:32:48.68 /ZlhhGjJ.net
>>555
レスありがとうございました。
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
578:１３２人目の素数さん
18/11/08 23:35:10.38 /ZlhhGjJ.net
>>557
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
579:１３２人目の素数さん
18/11/09 00:20:12.33 BZRFS9lT.net
不勉強ながらrle関数って初めて知ったけど使いやすそうだな
580:１３２人目の素数さん
18/11/09 00:36:46.75 Ui6BpICy.net
>>557
rleは0も数えているから間違っているんじゃ？
581:１３２人目の素数さん
18/11/09 02:01:36.95 Qla0VxTD.net
>>560
ご指摘ありがとうございました。
修正しました。
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
>
582:１３２人目の素数さん
18/11/09 07:25:00.97 Qla0VxTD.net
無理矢理１行にして実行
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
583:１３２人目の素数さん
18/11/09 07:59:47.47 rijPDFSR.net
>>562
意味もなくforループ回してた上に毎回sum使って真偽値を数値に変換してたけど
replicate使って最後に一回だけmean取ると2.066→1.886で1割短くなるのね
他人のコード読むのは勉強になる
584:１３２人目の素数さん
18/11/09 09:06:56.99 Ui6BpICy.net
>>561
whichいらねーよ。
585:１３２人目の素数さん
18/11/09 09:44:07.41 5AnUTlVm.net
>>561
全部が０のとき、エラーになるので修正
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
586:１３２人目の素数さん
18/11/09 10:10:01.39 Ui6BpICy.net
>>565
センスないなあ。
rle01 <- function(x) {
r <- rle(x)
one <- r$values == 1
if (any(one)) max(r$length[one]) else 0
}
587:１３２人目の素数さん
18/11/09 10:40:31.84 5AnUTlVm.net
>>566
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
>
588:１３２人目の素数さん
18/11/09 10:52:23.39 Ui6BpICy.net
>>567
そんな特別な場合のことに対して高速化するのは愚の骨頂。
1が一つでもある場合にsumを呼ぶのは余計じゃないのか？
余計なwhichがあるくらいだから、まずは素直にやりたいこと、やるべきことを正しく書くようにしたら？
589:１３２人目の素数さん
18/11/09 11:01:40.74 5AnUTlVm.net
1000本に1本あたる宝クジを100本買って続けて2本あたる確率のシミュレーション解の算出時間比較。
590:確率の理論値は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
591:１３２人目の素数さん
18/11/09 15:22:38.24 5AnUTlVm.net
>>569
分数で表すと
890788167367/9007199254740992
9.889735334744909e-05
592:１３２人目の素数さん
18/11/09 19:43:44.34 0M0agNOf.net
JPXのデータ、ファイル形式csvを読み込もうとするとうまく行かないんですが
どんな引数をつければいいですか
593:１３２人目の素数さん
18/11/10 10:36:49.41 QJ6NJqU7.net
>>568
実はシミュレーションじゃなくて
漸化式からのプログラム解を分数表示するプログラムはpythonで作成済。
ここに置いた。
URLﾘﾝｸ(egg.2ch.net)
594:１３２人目の素数さん
18/11/10 11:20:41.80 QJ6NJqU7.net
コインを100回ふったときの表連続の最大数が5であったときの
このコインの表がでる確率の期待値、モード比、信頼区間を求めるのが次のネタ。
unirootで算出できたけど
シミュレーションはどうすればいいのかアイデアが浮かばない。
MCMCで解決できるかなぁ？
これをシミュレーションで検証したい。
$Rscript main.r
lower mean mode upper
0.2487456 0.4469764 0.4589692 0.6386493
URLﾘﾝｸ(tpcg.io)
595:１３２人目の素数さん
18/11/11 20:11:21.81 ODPKEEGK.net
1億回のコイントスで何回連続して表がでる確率が高いかRでやってみた。
# 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万回で暴走。
ｽﾚﾘﾝｸ(hosp板:132番)
596:１３２人目の素数さん
18/11/12 21:04:57.19 boi8bvdM.net
"マラソン大会の選手に１から順番に番号の書かれたゼッケンをつける。
用意されたゼッケンN(=100)枚以下の参加であった。
無作為に抽出したM(=5)人のゼッケン番号の最大値はMmax(=60)であった。
参加人数推定値の期待値とその95%
597:信頼区間を求めよ" 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 結果は確認できたけど、もっと高速なシミュレーションアルゴリズムはあるだろうか？
598:１３２人目の素数さん
18/11/16 13:44:59.80 U19cHKqd.net
重複があるか否かを返す、anyDuplicatedという関数を知ったので総当たり比較と早いかどうか比べてみた。
覆面算を □/□ ＊ □/□ ＝ □□/□を解くの使ってみた。
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でコードは短くなるが速さが犠牲になった。
599:１３２人目の素数さん
18/11/16 16:13:49.03 zcEc4+wx.net
カルマンフィールタ
600:１３２人目の素数さん
18/11/18 12:56:58.01 tDkQkz0a.net
初歩的な質問で恐縮ですが、
truehistで出したヒストグラムの情報をテキストファイルとして保存する方法、教えて頂けませんか？
601:１３２人目の素数さん
18/11/18 16:15:11.61 hHA2tyUi.net
>>578
truehist関数はhist関数と違って戻り値を返さないから無理なのでは
602:１３２人目の素数さん
18/11/18 17:12:37.03 HQweHJGH.net
>>578
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"
で出力されるけど、
どのパラメータがあればヒストグラムが再現できるのかは不勉強に