19/12/14 04:10:22.87 L1XaepqW.net
シミュレーションその３ （処理速度の関係で１０万回の平均）
gacha=c(rep(1,36),rep(0,1000-36))
sim3 <- function() any(replicate(10,sample(gacha,1))==1)
mean(replicate(1e5,sim3()))
> mean(replicate(1e5,sim3()))
[1] 0.30691
これも0.307弱だな。
何が悪いんだろ？ 俺の頭かな？？
782:１３２人目の素数さん
19/12/14 05:45:07.83 qYB5MEbs.net
3.6%のガチャを10回回して全部外れる確率は
(1 - 0.036) ^ 10 ≒ 0.6930592
したがって、
3.6%のガチャを10回まわして1回以上当たる確率は
1 - 0.6930592 = 0.3069408
783:１３２人目の素数さん
19/12/14 08:18:13 L1XaepqW.net
>>747
ありがとう、
やっぱ、
これだったｗ
何が悪いんだろ？ 俺の頭かな？？
784:１３２人目の素数さん
19/12/14 08:19:41 L1XaepqW.net
>>748
シミュレーションにバグはなかったｗ
785:１３２人目の素数さん
19/12/14 12:29:26.79 cW38CJV+.net
小学生でも解けるぞ。がんばれ。
786:１３２人目の素数さん
19/12/14 16:29:58.85 QSWNkRNd.net
確率統計なんぞ無意味つまりインチキ
事実、結果こそすべて
787:１３２人目の素数さん
19/12/14 19:21:15.54 1VaGaO0p.net
その確率計算の激ナイーブな解法を示す
激レアは、レアだから1個とみなす。
故に、母集団の個数は
1÷0.036 = 27.777… きっと28個だ
手順1) 28枚のカードがある
手順2) 重複しない1~28の番号を振る
手順3) 28枚のカードをシャッフル
手順4) 1~10枚目のどれか1となる確率
と絶対同じハズ、だから、
P(1枚目で当) =27P27 ÷ 28P28 = 1/28
P(2枚目で当) も同様に1/28
P(3枚目で当) も同様に1/28
…
P(10枚目で当) も同様に1/28
で此等10個の事象は背反事象だから、
P = 10/28 = 0.357
∵有効数字3桁と勝手にしちゃう
ところてガチャってゲーム何か
よくわかんないけど、計算しちゃった
788:１３２人目の素数さん
19/12/14 20:44:12.02 1VaGaO0p.net
追記というか突然ですが、
そのガチャ3.6%の件、
超幾何分布なのか。二項分布なのか。
確率の小さいとか、母集団が小さい とかだと無視できないと思われる。
確率統計はギャンブル派生数学ぢゃ。
生半可な知識ではカモにされる。
現代の若者たちは、数学特に
確率統計はじめとするギャンブル
の能力が特段に欠けており、
R言語等のプログラミング教育で
ギャンブルゲームを学習すべきだ。
奇麗事の学問だけの今日の数学ぢゃ
カモにされるだけ。
健全な娯楽として賭博系確率統計学
をC R Java Pyson Javascript BASIC
の何れかを学校で学習すべきだ。
789:１３２人目の素数さん
19/12/16 09:18:31.40 t3hxiPBr.net
私のガチャのイメージは、
かつて、昭和の駄菓子屋によくある
ガチャガチャすなわちカプセルトイ
そのような健全的なギャンブルが
超幾何分布とか二項分布の理解に
役立つのだ。
限られたお小遣の10円玉数枚で、
如何にレアアイテムをGetするかを
子供らは、思考するからだ。
「残り物には福がある」は
確率統計的には正しいのか
子供同士で文学的に議論したものだ。
さて、今のガチャは恐らくは、
デジタルの媒体のスマホゲームだ。
二項分布でよいだろう。
R言語には、
二項分布の密度関数、それの累積関数
はモチロン、それに従う乱数生成を
提供してるようだ。
ゲームの仕組みが複雑化する今、
R言語等の乱数生成プログラミングで
これからのデジタル化ギャンブル社会
で、お金より大切な激レアをドンドン
無限にゲットできる人材の増加を
期待できる可能性を秘めている
790:１３２人目の素数さん
19/12/17 06:37:18.85 gnX7/8eN.net
カプセルトイ（ガチャ）１台に1000個カプセルが入っていて36個がアタリ(レアアイテム）とする。
同じカプセルトイが10台ある。カプセル取り出し後は補充されない。
アタリを１個でも手に入れる確率は
１台から10個取り出す場合(G10)と
１台から１個を１０台で取り出す場合(G01)
ではどちらが高いか？
そのシミュレーション
rm(list=ls())
N=1000 ; K=36 ;n=10 # アタリ3.6％
g=rep(c(1,0),c(K,N-K))
G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元（超幾何分布）
G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
mean(replicate(1e6,G10()))
mean(replicate(1e6,G01()))
791:１３２人目の素数さん
19/12/17 08:10:38.97 gnX7/8eN.net
シミュレーションと理論値
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.307745
[1] 0.3081121
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.307295
[1] 0.3069408
１台から１０個取り出した方がいいみたい。
792:１３２人目の素数さん
19/12/17 14:14:51.68 JPqkotiS.net
>>756 なるほど、
どの台も全部で、1000カプセルで、
どの台も当りが、36カプセルだと
G10(1台で10個)取り出す方が
僅かですが、確率良さそうですね。
シミュレーションも理論値も
同様な結論のようであり、
G10(1台で10個)取り出す方が僅かに
有利が分かり何か楽しかったです。
仮に、もしどの台も
250カプセル中当たり9カプセルなら、さらにG10戦略がG01戦略より有利な
感触を掴めました。
793:１３２人目の素数さん
19/12/17 15:44:13.73 gnX7/8eN.net
>>757
1000個だとシミュレーションの差が微妙で再現性が不安だったけど
250個にすると、差が明らかにつきますね。
> N=250 ; K=9 ;n=10 # アタリ3.6％
> g=rep(c(1,0),c(K,N-K))
>
> G10 <- function() any(sample(g,n,replace=FALSE)>0) # 非復元（超幾何分布）
> G01 <- function() any(sample(g,n,replace=TRUE )>0) # 復元(二項分布)
>
> mean(replicate(1e6,G10())) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.31082
[1] 0.3117069
> mean(replicate(1e6,G01())) ; 1-(1-K/N)^n
[1] 0.306782
[1] 0.3069408
794:１３２人目の素数さん
19/12/18 14:18:29.16 GQllmBub.net
rhyper で超幾何分布の乱数を発生させることができたんだな。
これと rbinomを使うと sample関数を使わなくても同じことができた。
N=250 ; K=9 ;n=10 # アタリ3.6％
mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n) # 非復元（超幾何分布）
mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n # 復元(二項分布)
結果
> mean(replicate(1e6,rhyper(1,K,N-K,n)>0)) ; 1-choose(N-K,n)/choose(N,n)
[1] 0.311581
[1] 0.3117069
> mean(replicate(1e6,any(rbinom(n,1,K/N)>0))) ; 1-(1-K/N)^n
[1] 0.307213
[1] 0.3069408
1台から順次１０個取り出すのをイメージすれば残り物には福があるとも言えなくもないな。
>
795:１３２人目の素数さん
19/12/18 21:18:37 GQllmBub.net
"フランスの数学者パスカル（1623～1662)が1654年にフェルマーにあてた手紙が、現在の確率
論の始まりだと言われている。当時の有名な賭博師メレがパスカルに以下のような問題を持ち
込み、その問題についてがその手紙のやりとりの中で論じられているそうである。
甲乙二人がおのおの32ピストル（当時のお金の単位）の金を賭けて勝負したとする。
そしてどちらかが先に３点を得たものを勝ちとし、勝った方がかけ金の総額64ピストルをもら
えるとする。ところが甲が２点、乙が１点を得たとき、勝負が中止になってしまった。
このとき、二人のかけ金の総額64ピストルを甲と乙にどのように分配すればよいだろうか。
ただし二人の力は互角で、勝つ確率はそれぞれ1/2ずつだとする。"
# 先にw(=3)点を得たものを勝ち、甲がA(=2)点、乙がB(=1)点を得たとき、勝負が中止
# 甲の確率を求める
gambling <- function(A=2,B=1,w=3,k=1e5){ # k : how many simulate
sim <- function(){
while(A < w & B < w){
g = rbinom(1,1,p=0.5)
if(g==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}
> gambling(2,1,3)*64 #64ピストル配分
[1] 48.13056
> gambling(2,0,4) # 日本シリーズ
[1] 0.80945
796:１３２人目の素数さん
19/12/19 09:59:55.45 V+OT4hGF.net
>>760の計算結果は、興味深い結果だ
で、
自己流の直感でラフな計算をしてみた。
甲1点、乙1点 ⇒1:1按分 ∵自明 ★
甲3点、乙1点 ⇒1:0按分 ∵自明 ☆
安直に★と☆の平均を考えると
甲2点、乙1点 ⇒1:0.5按分 ∵安直
で、総額で64ピストルだから、
甲に、64*1/1.5 = 42.67 ピストル
乙に、64*0.5/1.5 = 21.33 ピストル
を配分すれば、思ってしまった。
ラフな計算だと、少々不公平なことが
起こるようだ。
R言語等のシミュレーションは、
超天才のワィの直感より、
さらに厳密な確率計算ができそうだ。
797:１３２人目の素数さん
19/12/19 16:21:08.12 kVFzAP55.net
>>760
最後はAが得点して勝負が決定するのを忘れずにAが勝者になるまでのBの点数で場合分け
NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B：現在の点数 ,p:甲の勝率
ans=0
for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数
ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A)
}
return(ans)
}
> NS(3,2,1,0.5)*64 #64ピストルの分配
[1] 48
> NS(4,1,0,0.5) #日本シリーズで先勝したほうが優勝する確率
[1] 0.65625
NSは日本シリーズの頭文字から命名、No Skinではないｗｗ
798:１３２人目の素数さん
19/12/19 18:00:34.70 V+OT4hGF.net
賭博等の確率計算をイメージすると
ワィの脳内は
799:活性化される だから、 ドンドン、妄想全開となるぅぅぅ。 という訳で、妄想内容を記載する 17世紀の地球に行き、胴元になって 賭博ビジネスでガッポリ儲けるゼェ。 17世紀の地球で、胴元になって、 コイントス賭博事業をするゼェ。 で、システムは、2回コイントスし、 表の回数が0回⇒参加者が32万円GET 表の回数が1回⇒胴元に、32万円PAY 表の回数が2回⇒参加者が32万円GET 参加料は、1万円 まんまと、参加者が沢山 1000人いた ある参加者は、以下の様に怪説をした 組み合わせは、3通りだ。 1通り目 表の回数が0回 2通り目 表の回数が1回 3通り目 表の回数が2回 だから、参加者の勝つ確率は2/3だ で、モチロン、胴元ワィの商売は成功 およそ、1000万円儲かっちゃた。
800:１３２人目の素数さん
19/12/20 02:07:38.50 yiLw1Jz8.net
0745
しろ@huwa_cororon 11月27日
苦節6ヶ月、初満点&一等賞です！
URLﾘﾝｸ(twitter.com)
(deleted an unsolicited ad)
801:１３２人目の素数さん
19/12/21 10:49:06.51 +RV2yvS7.net
日本シリーズで賭けをする。
日本シリーズは先に4勝したチームが優勝。
勝率はそれまでの（引き分けを除いた）通算勝率に従うとする。
シリーズ開始前の通算成績はA:2勝、B:4勝であった。
今シリーズでAが先勝（第一試合に勝利）した。
この時点でA,Bのどちらに賭ける方が有利か？
802:１３２人目の素数さん
19/12/21 22:03:14.56 vLicX+v2.net
分からない問題スレで３進法の小数とかが話題になっていた
ｽﾚﾘﾝｸ(math板:845番)
>>
845 名前：１３２人目の素数さん[] 投稿日：2019/12/21(土) 05:52:50.74 ID:wum1jR1j
...
10進数で表された5/9を3進数になおせという問題があって、答は0.12だったのですが、
...
<<
こんなのを作ってみた。
# 小数点付きの数numをN進法で表示する
rm(list=ls())
dec2basen <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納
x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
if(N<=10){ # Nが10以下は数値として表示
cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='')
}
else{ # Nが1１以上は整数部分と小数部分を数列で表示
print(list(int=r,decimal=a))
}
invisible()
}
> dec2basen(5/9,3)
0.120000000000000
> dec2basen(3.14159265359,7)
3.06636514320
> dec2basen(3.14159265359,8)
3.11037552421
> dec2basen(3.14159265359,9)
3.12418812407
> dec2basen(3.14159265359,14)
$int
[1] 3
$decimal
[1] 1 13 10 7 5 12 13 10 8 1 4
803:１３２人目の素数さん
19/12/22 08:18:39.30 Ez3wWboE.net
>>766
バグの原因を追求したら、Rの仕様のせいみたいだな。
> (1.2-1)*5==1
[1] FALSE
> (1.2-1)*5>1
[1] FALSE
> (1.2-1)*5<1
[1] TRUE
これがfloor関数が誤値を返してくる原因だった。
> floor(1)
[1] 1
> floor((1.2-1)*5)
[1] 0
804:１３２人目の素数さん
19/12/22 08:52:46.43 UYCjcQMO.net
>>767
数値演算についてもっと勉強しな。
たいていの言語でそうなる。
805:１３２人目の素数さん
19/12/22 13:06:36.02 Ez3wWboE.net
すると、pythonに移植しても同じ結果ということか。
IPython 6.5.0 -- An enhanced Interactive Python.
C:\bin\Anaconda3\lib\site-packages\ipykernel\parentpoller.py:116: UserWarning: Parent poll failed. If the frontend dies,
the kernel may be left running. Please let us know
about your system (bitness, Python, etc.) at
ipython-dev@scipy.org
ipython-dev@scipy.org""")
(1.2-1)*5==1
Out[1]: False
(1.2-1)*5>1
Out[2]: False
(1.2-1)*5<1
Out[3]: True
(1.2-1)*5
Out[4]: 0.9999999999999998
806:１３２人目の素数さん
19/12/22 17:15:43.69 Ez3wWboE.net
round を使って回避することにした
> 0.728*5-3.64
[1] -4.440892e-16
> round(0.728*5-3.64,5)
[1] 0
んで、デバッグ版
# 小数点付きの数numをN進法で表示する(62進法まで対応 0-9,a-z,A-Z)
rm(list=ls())
dec62 <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示数列を格納
k=max(nchar(num)-nchar(floor(num))-1,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
x=round(num-floor(num),k) # e.g. 7.28-floor(7.28)-0.28 != 0に対応
for(i in 1:k){
y=round(x*N,k) # e.g. 0.728*5-3.64 !=0 に対応
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters,LETTERS)[1:N]
if(N<=62){ # Nが62以下は数値として表示
cat(paste(fig[b$integer+1],sep='',collapse=''),
'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが63以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b) # b$num:検証用
}
807:１３２人目の素数さん
19/12/22 17:25:01.53 Ez3wWboE.net
>>769
Haskellでもやってみた
Prelude> (1.2-1)*5
0.9999999999999998
Prelude> 0.72*5-3.6
-4.440892098500626e-16
808:１３２人目の素数さん
19/12/25 07:12:30.44 oEKznZ6+.net
>>768
レスありがとうございました。
アドバイスに従ってちょっと勉強してみました。
他の言語に移植しても無駄とわかって助かりました。
こうなる理由が理解できました。
(1+1/10-1)*10==1
[1] FALSE
> (1-1+1/10)*10==1
[1] TRUE
809:１３２人目の素数さん
20/01/06 12:56:45.01 e9wyGMBv.net
からnまでの順列を列挙するスクリプト。
これをwhile使って高速化するにはどうすればいいだろう？
俺の頭では思いつかない。
perm <- function (n) {
v=1:n
sub <- function(n, v) {
if (n == 1)
matrix(v, 1, 1)
else {
x = NULL
for (i in 1:n) x =rbind(x, cbind(v[i], sub(n -1,v[-i])))
x
}
}
sub(n, v)
}
> perm(4)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 1 2 4 3
[3,] 1 3 2 4
[4,] 1 3 4 2
[5,] 1 4 2 3
[6,] 1 4 3 2
[7,] 2 1 3 4
[8,] 2 1 4 3
[9,] 2 3 1 4
[10,] 2 3 4 1
[11,] 2 4 1 3
[12,] 2 4 3 1
[13,] 3 1 2 4
[14,] 3 1 4 2
[15,] 3 2 1 4
[16,] 3 2 4 1
[17,] 3 4 1 2
[18,] 3 4 2 1
[19,] 4 1 2 3
[20,] 4 1 3 2
[21,] 4 2 1 3
[22,] 4 2 3 1
[23,] 4 3 1 2
[24,] 4 3 2 1
>
810:１３２人目の素数さん
20/01/06 21:34:47.52 yym51Tg7.net
>>773
自己解決 e1071というパッケージに再帰呼び出しなしのpermutationがありました。
library(e1071)
permu <- function (n) {
z <- matrix(1)
for (i in 2:n) {
x <- cbind(z, i)
a <- c(1:i, 1:(i - 1))
z <- matrix(0, ncol = ncol(x), nrow = i * nrow(x))
z[1:nrow(x), ] <- x
for (j in 2:i - 1) {
z[j * nrow(x) + 1:nrow(x), ] <- x[, a[1:i + j]]
}
}
return(z)
}
811:１３２人目の素数さん
20/01/09 14:28:06.78 rRnLkyt7.net
可変引数の扱い
Cで
#include <stdio.h>
int main(int argc, char* argv[])
{
while(argc--)
printf("%s\n", *argv++);
return 0;
}
をRでやるには ... を使って、c(...)もしくはlist(...)で引き出せばいいんだな。
検索してもなかなかわからなかった。
main <- function(...){
argv=c(...)
for(i in 1:length(argv)) cat(argv[i],'\n')
}
812:１３２人目の素数さん
20/01/15 18:08:54 l229ykjv.net
>>548
library(Rmpfr)
mpfr(factorial(50),5000)
one = mpfr(1, 5000)
factorial(one*50)
> mpfr(factorial(50),5000)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713018969967457666435945132957882063457991132016803840
> one = mpfr(1, 5000)
> factorial(one*50)
1 'mpfr' number of precision 5000 bits
[1] 30414093201713378043612608166064768844377641568960512000000000000
何故か、oneを高精度に設定しておくと、正しい計算をしてくれた。
813:１３２人目の素数さん
20/01/26 16:33:21 G7gVG9Ku.net
RのroundのIEEE仕様
> round(2.5)
[1] 2
> round(1.5)
[1] 2
は丸め誤差が減るという説明だったので実験してみた。
# 四捨五入 vs ＩＥＥＥ with FUN(=mean) for x(=c(-0.9,-0.8,...,0.8,0.9))
comp <- function(n=10,FUN=mean,x=sample((-9:9)/10,n,rep=T),print=T){
X=FUN(x)
dif = abs(X-FUN(round(x))) - abs(X-FUN(f45(x))) # round後の実行と四捨五入後の実行の差
if(dif!=0 & print) cat(dif<0,' : ',sort(x),'\n') # 差があれば表示 round後が小さければTRUE
return(dif) # 差を返す dif<0:round優位 dif>0:四捨五入優位
}
comp()
k=1e5
# mean
re=replicate(k,comp(print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
# squared sum
comp(FUN=function(x) sum(x^2))
re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
> # 平均 mean
> re=replicate(k,comp(print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ 四捨五入 引き分け
0.25235 0.16199 0.58566
> # 平方和 squared sum
> re=replicate(k,comp(FUN=function(x)sqrt(sum(x^2)),print=F))
> c(ＩＥＥＥ=mean(re<0),四捨五入=mean(re>0),引き分け=mean(re==0))
ＩＥＥＥ 四捨五入 引き分け
0.40519 0.01413 0.58068
平均だと小差だが、統計計算の内部処理で使う平方和だと大差がついたので、体感できた。
814:１３２人目の素数さん
20/02/08 17:42:30.00 +JttXwIS.net
"隔離中のクルーズ船では船内の換気が共通らしいから運悪く13日後に発症した奴がいるとその近くの部屋のやつは
プラスで14日隔離しないといけない。それが今の船内の状況という。
両隣のどちらかが感染したら14日延長、どの部屋も1日で感染する確率pは1%
部屋の配置は長方形でどの部屋にも両隣があるとする。
発症するか、隔離期間が終われば下船できる。
全員定員１の個室として客と乗務員を合わせた人数nは3000人。
クルーズ船から全員下船できる日数の期待値は？"
815:１３２人目の素数さん
20/02/08 17:42:46.53 +JttXwIS.net
rm(list=ls())
p=0.001 # 1日の感染確率
n=3000 # 収容者人数
d14=14 # 隔離日数
Q0=rep(d14,n) # 必要隔離日数の配列の初期値
#
ODA <- function(X){ # One Day After 必要隔離日数の配列 x から１日後の配列を返す
Q=X #
816:Xを作業用Qに代入 S=which(Q!=0) # 未感染者susceptibleのindex s=length(S) # 未感染者人数 m=rbinom(1,s,p) # その日の感染数 if(m==0){ # 感染者0なら未感染者の Q[S]=Q[S]-1 # 必要隔離日数を１日減らして返す invisible(Q) }else{ I.idx=sample(S,m) # 未感染者Sからm人が感染、感染者のindex (Infected.idx) id2ext <- function(id){ # 感染者idからの両隣のidを返す plus =ifelse(id==n,1,id+1) # 最大番号なら1を返す,それ以外は1増やす minus=ifelse(id==1,n,id-1) # １なら最大番号を返す,それ以外は１減らす unique(c(plus,minus)) # 重複を除く } E.idx=as.vector(sapply(I.idx,id2ext)) # 感染者の両隣のidを列挙し E.idx=unique(E.idx) # 重複を除く Extended quarantine.idx Q[E.idx]=d14 # 両隣は必要隔離日数をd14にリセット Q[I.idx]=0 # 感染者は必要隔離日数=0 (リセット者に重複していたら０で上書き) IE.idx=unique(c(I.idx,E.idx)) # 感染もしくは感染者の隣のindex minus1ifnotzero <- function(x) ifelse(x==0,0,x-1) # ゼロでなければ１を減じて返す Q[-IE.idx]=minus1ifnotzero(Q[-IE.idx]) # 感染もなく感染者の隣りでもない人の隔離日数を１日減らす invisible(Q) } } sim <- function(){ AZ=FALSE # All Zero 隔離日数がすべて０か？ i=0 # カウンター Q=Q0 # 隔離日数配列初期値 while(!AZ){ # 隔離日数がすべて0になるまで Q=ODA(Q) # One Day After関数を実行 AZ=all(Q==0) # All zero？ i=i+1 # カウンター } return(i) # 何日かかったかを返す } k=1e4 # k回のシミュレーションの RE=replicate(k,sim()) mean(RE) # 平均値＝期待値
817:１３２人目の素数さん
20/02/08 18:25:56.61 Xaabj2Ml.net
感染確率0.1％やん
低すぎぃ
818:１３２人目の素数さん
20/02/08 19:50:08.59 +JttXwIS.net
>>780
一日の感染確率。
3000人いて1日3人くらい報告されたから0.1%/日くらいかな
819:１３２人目の素数さん
20/02/08 19:51:16.44 +JttXwIS.net
実行したら全員が下船できるまで約1ヶ月になった。
コードにまだバグがあるかも。
820:１３２人目の素数さん
20/02/09 21:31:53 IDQ3+Iu7.net
Rのcsvの文字化け治らねえ
csvそのものをutf-8にしたら今度は不正なマルチバイトがあるってエラー
どうすればいいか教えて
821:１３２人目の素数さん
20/02/09 22:07:23 vcjg+5qP.net
>>783
どんなデータ？
ヘッダーだけじゃなくてセル内も日本語？
822:１３２人目の素数さん
20/02/09 23:10:12 KodPd2ez.net
>>783
エスパーじゃないから、まず確認だけど
・Rを使ってる環境(OS)とRのバージョンは？
・CSVのコード変換に使ってるソフトは？
・CSVを読み込んでる関数と指定している引数は？
823:１３２人目の素数さん
20/03/06 19:53:26 WwE0ymvi.net
nCov2019 for studying COVID-19 coronavirus outbreak URLﾘﾝｸ(guangchuangyu.github.io)
後で遊んでみよう
824:１３２人目の素数さん
20/03/07 09:46:33.43 3M7vmA2q.net
>>786
エラーがでた。R 3.6.3にアップデートしたけど
> library('remotes')
> remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
Downloading GitHub repo GuangchuangYu/nCov2019@master
Skipping 1 packages not available: BiocStyle
Installing 5 packages: downloader, cowplot, maps, magick, BiocStyle
Error: Failed to install 'nCov2019' from GitHub:
(converted from warning) package ‘BiocStyle’
825:is not available (for R version 3.6.3)
826:１３２人目の素数さん
20/03/07 11:22:48.00 tNSvsFaR.net
>>787
エラーと通り、BiocStyleパッケージがないので別途、bioconductorからインストールすればいい。
URLﾘﾝｸ(bioconductor.org)
827:１３２人目の素数さん
20/03/07 11:34:31.80 3M7vmA2q.net
>>788
ありがとうございます。
828:１３２人目の素数さん
20/03/07 11:38:15.72 3M7vmA2q.net
We also developed a web app (URLﾘﾝｸ(www.bcloud.org)) with interactive plots and simple time-series forecasts.
という記述があるんだけど 関数 nls による非線形回帰分析 で作成しているのかなぁ？
829:１３２人目の素数さん
20/03/09 20:05:18 0N1NTePA.net
面白い問題スレをRでシラミ潰しにカウントしてみた。
"縦n個、横n個のマス目のそれぞれに 1,2,3,...,n の数字を入れていく。
このマス目の横の並びを行といい、縦の並びを列という。
どの行にも、どの列にも、2つの対角線上にも同じ数字が1回しか現れない入れ方は何通りあるか求めよ。(2020京大文系 改)
n=4がオリジナル
"
library(gtools)
n=4
(a=permutations(n,n)) # ｎの順列
r=nrow(a)
f<-function(x){ # x=c(1,2,3,4) -> rbind(a[1],a[2],a[3],a[4])
ans=NULL
for(i in x){
ans=rbind(ans,a[i,])
}
return(ans)
}
b=permutations(r,n)
head(b) ; tail(b)
B=apply(b,1,f)
diag2 <- function(x){ # 行列xの右上から左下への対角線の配列を返す
n=nrow(x)
ans=numeric(n)
for(i in 1:n) ans[i]=x[i,n+1-i]
return(ans)
}
g <- function(v){ # matri(v,n)で列、対角線の要素がすべて異なればTRUEを返す
n=sqrt(length(v))
(x=matrix(v,n))
if(length(unique(diag (x))) < n) return(FALSE)
if(length(unique(diag2(x))) < n) return(FALSE)
flag=TRUE
for(i in 1:n){
if(length(unique(x[,i])) < n){
flag=FALSE
break
}
}
return(flag)
}
counter=NULL
for(i in 1:ncol(B)){
if(g(B[,i])) counter=append(counter,i)
}
length(counter)
counter
matrix(B[,counter[1]],n)
matrix(B[,counter[2]],n)
matrix(B[,counter[48]],n)
830:１３２人目の素数さん
20/03/10 07:11:26.36 H1fx2jVB.net
>>791
n=5 にするとメモリー不足でエラーがでた。
サンプリングでやってみる
# simulation
n=5
sim <- function(n=5){
v=as.vector(t(replicate(n,sample(n)))) # n × n行列要素をベクトル化
if(g(v,dia=TRUE)){ # diagonal latin squareなら
print(matrix(v,n)) # 行列化して表示
return(1)
}
return(0)
}
k=1e6
sum(replicate(k,sim(5)))
831:１３２人目の素数さん
20/03/10 12:07:20.32 KHQs8ybj.net
１億回モンテカルロしてようやく、1個みつかった。
[,1] [,2] [,3] [,4] [,5]
[1,] 2 5 4 3 1
[2,] 4 3 1 2 5
[3,] 1 2 5 4 3
[4,] 5 4 3 1 2
[5,] 3 1 2 5 4
832:１３２人目の素数さん
20/03/11 11:06:12 7KW6Eui5.net
嫌ね最近は総当たりアルゴリズムで枝刈り使えば実用的な時間で解決できるこ知らない人が多くて
833:１３２人目の素数さん
20/03/15 08:45:11.95 OTl1KJku.net
versionが3.6,xになって、str2lang とか str2expressionという関数が追加されたんだな。
URLﾘﾝｸ(id.fnshr.info)
で知った
例
rep("1:6",5)
paste(rep("1:6",5),collapse=',')
str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))
eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')))
(s=rep("1:6",5))
(str=paste(s,collapse=','))
(lang=str2lang(paste("expand.grid(",str,')')))
eval(lang)
834:１３２人目の素数さん
20/03/21 17:23:30 RyI2Q/uv.net
The Clinical and Chest CT Features Associated with Severe and Critical COVID-19 Pneumonia
URLﾘﾝｸ(journals.lww.com)
の付表にこんなデータがあった
# Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value
# Number of lobes Median (interquartile range, IQR)
# 5 (4-5) 5 (5-5) 5 (3-5) 0.003
重症群２５例 中等症群５８例でCTで肺炎像が確認できた肺葉数(１～５）が
全体で中央値が５ 第一四分位数(少ない方から数えて25%の値）が４、第三四分位数(少ない方から数えて75%の値）が５
重症群では中央値５、第一四分位数５、第三四分位数５
中等症群では中央値５、第一四分位数３、第三四分位数５で
p値が０．００３であるという。
０．００２５～０．００３５までとして、
この条件に見合う８３例をリストアップすると
> re[100,][1:25]
[1] 1 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
> re[100,][26:83]
[1] 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[48] 5 5 5 5 5 5 5 5 5 5 5
という類の組み合わせで
19504通りと計算された。
重症群の肺葉数の平均は４．７６
中等症群の平均は４．０４
と計算された。
タイ値がこれだけ多いデータにWilcox検定を使う意味がわからんな。
835:１３２人目の素数さん
20/03/22 08:57:06 liILqu/N.net
# 感度0.7 特異度0.9のPCR検査で1000人中10人が陽性だったときの有病率の信頼区間をシミュレーションして計算するスクリプト
# アルゴリズムには自信がないｗ
sim <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){
# n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画
prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定
PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算
prev.post = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算
# NPVでも同じ結果
# NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算
# prev.post = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算
mean=mean(prev.post) # 期待値
median=median(prev.post) # 中央値
mode=density(prev.post)$x[which.max(density(prev.post)$y)] # 最頻値
LU=unlist(HDInterval::hdi(prev.post))[1:2] # 95%信頼区間
re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値
if(print){ # ヒストグラム描画
par(mfrow=c(2,2))
hist(prev,freq=F,xlim=c(0,1)) ; lines(density(prev),lwd=2)
hist(PPV,freq=F,xlim=c(0,1)) ; lines(density(PPV),lwd=2)
hist(prev.post,freq=F,xlim=c(0,1)) ; lines(density(prev.post),lwd=2)
BEST::plotPost(prev.post)
par(mfrow=c(1,1))
}
print(round(re,4)) # 概算値表示
invisible(list(re,prev.post)) # 代表値と事後有病率を非表示で返す
}
sim(100,1)
sim(1000,10)
sim(10000,100)
sim(100000,1000)
836:１３２人目の素数さん
20/03/22 14:13:54 liILqu/N.net
>>797(自己レス)
間違いに気づいたので撤回します。
837:１３２人目の素数さん
20/03/23 14:55:08 iIfIhG5+.net
MCMCで解決してみた。
【富山県最強伝説】新型
838:コロナウイルスPCR検査件数 54人 陽性0人 ある集団から54人を無作為に選んでPCR検査したら陽性0であった。 PCR検査の感度0.7 特異度0.9としてこの集団の有病率の期待値と95%信頼区間を推測せよ。 #---------------------------------------------------------------------------- # 感度SEN, 特異度SPCの検査でN人中X人が陽性であったときの推定有病率prevalence # 弱情報事前分布はprevalence ~ dunif(0,UL) , UL:上限 #---------------------------------------------------------------------------- PCRj <- function(N,X,UL=1,SEN=0.7,SPC=0.9,verbose=TRUE){ # UL:upper limit of dunif(0,UL) library(rjags) modelstring=paste0(' model { x ~ dbin(p,n) p <- prev*sen + (1-prev)*(1-spc) prev ~ dunif(0,',UL,') } ') if(verbose & UL!=1) cat(modelstring) writeLines(modelstring,'TEMPmodel.txt') dataList=list(n=N,x=X,sen=SEN,spc=SPC) jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList,quiet=TRUE) update(jagsModel) codaSamples = coda.samples( jagsModel , variable=c("prev","p"), n.iter=1e6, thin=10) js=as.matrix(codaSamples) BEST::plotPost(js[,'prev'],xlab='prevalence',showMode = TRUE) ; lines(density(js[,'prev']),col='skyblue') round(c(mean=mean(js[,'prev']),HDInterval::hdi(js[,'prev'])[1:2]),10) } PCRj(100,1) PCRj(1000,10) PCRj(10000,100) PCRj(100000,1000) PCRj(54,0) PCRj(54,0,UL=0.1)
839:１３２人目の素数さん
20/04/03 19:48:06 MNrUuJMd.net
a b c d e f g h i j
1 -0.0377710 0.02994835 0.1779426 0.29057472 0.13252796 0.3279709 -0.3534731 -0.54812874 -0.3869768 -0.1233967
2 0.2707033 0.09238235 -0.1042958 -0.08485129 -0.02590955 0.2159095 -0.3338677 0.00865952 -0.1575856 0.0953795
840:１３２人目の素数さん
20/04/03 19:50:07 MNrUuJMd.net
上のDFを
library(psych)
ICC(DF）
で級内相間求めようとすると変な値でるんだけどなんででしょう？
841:１３２人目の素数さん
20/04/03 22:01:48 1EdxzcOL.net
> DF
V1 V2
a -0.03777100 0.27070330
b 0.02994835 0.09238235
c 0.17794260 -0.10429580
d 0.29057472 -0.08485129
e 0.13252796 -0.02590955
f 0.32797090 0.21590950
g -0.35347310 -0.33386770
h -0.54812874 0.00865952
i -0.38697680 -0.15758560
j -0.12339670 0.09537950
> library(psych)
> ICC(x=DF,lmer=FALSE)
Call: ICC(x = DF, lmer = FALSE)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.36 2.1 9 10 0.13 -0.18 0.74
Single_random_raters ICC2 0.34 2.0 9 9 0.17 -0.26 0.74
Single_fixed_raters ICC3 0.32 2.0 9 9 0.17 -0.24 0.72
Average_raters_absolute ICC1k 0.53 2.1 9 10 0.13 -0.43 0.85
Average_random_raters ICC2k 0.51 2.0 9 9 0.17 -0.69 0.85
Average_fixed_raters ICC3k 0.49 2.0 9 9 0.17 -0.63 0.84
Number of subjects = 10 Number of Judges = 2
842:１３２人目の素数さん
20/04/03 23:40:27 MNrUuJMd.net
>>802
すいません
書き方が良くなかったです。
subject が１，２
judge a,b,c,dをきちんとfirst,second,third,fourth
と書くべきでした。
subjectとjudgeはあってると思うのです。
843:１３２人目の素数さん
20/04/03 23:42:39 MNrUuJMd.net
例えば
first second third fourth
0 0.05 0.1 0.15
0 0 0.05 0.05
0.15 0.1 0.05 0.1
0 0.05 0.2 0.05
0.1 0.15 0.1 0.05
0.15 0.05 0.1 0.05
というデータで
Call: ICC(x = over_sight)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.0068 1 5 18
844:0.43 -0.19 0.48 Single_random_raters ICC2 0.0068 1 5 15 0.44 -0.19 0.48 Single_fixed_raters ICC3 0.0068 1 5 15 0.44 -0.19 0.48 Average_raters_absolute ICC1k 0.0268 1 5 18 0.43 -1.70 0.79 Average_random_raters ICC2k 0.0268 1 5 15 0.44 -1.82 0.79 Average_fixed_raters ICC3k 0.0268 1 5 15 0.44 -1.82 0.79 Number of subjects = 6 Number of Judges = 4
845:１３２人目の素数さん
20/04/03 23:44:30 MNrUuJMd.net
となるのですが、ICC1が0.0068ってこんなに低くなるのでしょうか？
結果のとり得る値が1～0なので最大15％ぐらいの変動幅で、
制度としてはかなりいい検査方法だとおもったのですが・・・
846:１３２人目の素数さん
20/04/04 00:07:42 npCU99bV.net
少し理解できました。
範囲制約性の問題があるのですね。
今回のサンプルがみんな低い得点で分散が低いからICCが低く算出されていました。
ダミーデータ
1,1,1,1
を追加したらICC=0.99となりました。
得点の高いサンプルを追加しないとICCで信頼性は測れないのですね。
847:１３２人目の素数さん
20/04/04 00:09:17 npCU99bV.net
いま有病者と健常者の両方を計測する検査法で健常者のみのデータでICCを算出しようとしているのが間違いということですね。
上記のような場合で、健常者のみでも信頼性を算出する方法っていうのはないのでしょうか？
848:１３２人目の素数さん
20/04/04 21:44:33 ZFu90Xbq.net
>>806
DF=rbind(
c(0,0.05,0.1,0.15),
c(0,0,0.05,0.05),
c(0.15,0.1,0.05,0.1),
c(0,0.05,0.2,0.05),
c(0.1,0.15,0.1,0.05),
c(0.15,0.05,0.1,0.05),
c(0.0001,0.0001,0.0001,0.0001))
とダミーに0.0001を加えても
> DF2ICC1(DF)
ICC(1,1) ICC(1,k)
0.2462707 0.5665262
> psych::ICC(DF)
type ICC F df1 df2 p lower bound
Single_raters_absolute ICC1 0.25 2.3 6 21 0.072 -0.027
Single_random_raters ICC2 0.25 2.3 6 18 0.079 -0.028
Single_fixed_raters ICC3 0.25 2.3 6 18 0.079 -0.034
Average_raters_absolute ICC1k 0.57 2.3 6 21 0.072 -0.115
Average_random_raters ICC2k 0.57 2.3 6 18 0.079 -0.122
Average_fixed_raters ICC3k 0.57 2.3 6 18 0.079 -0.154
となるから、高値データでなくてもいいみたい。
849:１３２人目の素数さん
20/04/24 21:32:48 1Teq2ITX.net
R 4.0.0 来た
850:１３２人目の素数さん
20/04/25 06:37:54.89 R/UD6QQG.net
SIGNIFICANT USER-VISIBLE CHANGES
URLﾘﾝｸ(cran.r-project.org)
851:１３２人目の素数さん
20/04/25 07:20:54 R/UD6QQG.net
パッケージを一括インストール
install.packages(available.packages()[,1])
852:１３２人目の素数さん
20/04/25 07:29:33 R/UD6QQG.net
>>811
これやってからでないとパッケージのコンパイルができない。
After installation is complete, you need to perform one more step to be able to compile R packages:
you need to put the location of the Rtools make utilities (bash, make, etc) on the PATH.
The easiest way to do so is create a text file .Renviron in your Documents folder which contains the following line:
PATH="${RTOOLS40_HOME}\usr\bin;${PATH}"
You can do this with a text editor, or you can even do it from R like so:
writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")
853:１３２人目の素数さん
20/04/25 08:56:31 R/UD6QQG.net
>>811
libraryのフォルダで dir /B とやって既にinstall済のリストを作って
install.packages("rstan","ggplot2",....)とやった方がいいな。
5000ほどパッケージがあるらしいけど、全部インストールしたらどれくらいの容量が必要なんだろう？
854:１３２人目の素数さん
20/04/25 12:31:15 uVYqRSBX.net
やってみてくれｗ
855:１３２人目の素数さん
20/04/26 21:05:35 tZeixXMR.net
>>813 これでいいかな？
targetPackages <- available.packages()[,"Package"]
newPackages <- targetPackages[!(targetPackages %in% installe
856:d.packages()[,"Package"])] install.packages(newPackages, repos = "http://cran.r-project.org")
857:１３２人目の素数さん
20/04/26 22:22:15 tZeixXMR.net
高血圧で不安だという老人の経過観察入院1件。発熱してなくてほっとした。
ステーキハウスがテイクアウトを始めたからインセンティブはその足しにしようっと。
858:１３２人目の素数さん
20/05/03 15:15:43 ICXX0aJG.net
>>814
コンパイルなしでやってみた。
サイズ：19.6 GB (21,082,530,395 バイト)
ディスク上のサイズ：20.3 GB (21,805,858,816 バイト)
ファイル数: 527,982、フォルダー数: 118,521
859:１３２人目の素数さん
20/05/03 15:37:44 f6MQMkc5.net
Ubuntu 20.04でRcmdrがビルドできないんですけど、これはソフトウェアが20.04に対応するまで待つしかないですか？
860:１３２人目の素数さん
20/05/03 20:21:16 l1OdZntk.net
>>818
自己解決しました。ライブラリの不足でした。
861:１３２人目の素数さん
20/05/04 17:06:14.69 HJhKUMfp.net
4.0でまたsjis関連の問題が出まくりそう
862:１３２人目の素数さん
20/05/07 01:08:38 VnQvkZ57.net
4.0になってから
col=番号指定での色が目に優しくなったみたい。
hist(runif(100),col=2)
hist(runif(100),col=rgb(1,0,0,1))
hist(runif(100),col='red')
863:１３２人目の素数さん
20/05/07 18:12:52.58 nLAcI1eZ.net
rgbで256色でそんなに変化ないんじゃない。
864:１３２人目の素数さん
20/05/07 19:16:11 VnQvkZ57.net
>>822
colours() と入力すると 657色から選べるから256以上ありそう。
865:１３２人目の素数さん
20/05/07 23:34:33 8IDmAqjy.net
4.0からパレット変えたんじゃなかったか？
866:１３２人目の素数さん
20/05/07 23:39:55 8IDmAqjy.net
リリースノートより
The palette() function has a new default set of colours (which are less saturated and have better accessibility properties).
There are also some new built-in palettes, which are listed by the new palette.pals() function.
These include the old default palette under the name "R3". Finally, the new palette.colors() function allows a subset of colours to be selected from any of the built-in palettes.
867:１３２人目の素数さん
20/05/07 23:46:21 8IDmAqjy.net
具体的な変更点はこのあたりを参照方
URLﾘﾝｸ(developer.r-project.org)
868:１３２人目の素数さん
20/05/08 00:23:14.39 CfFk/Uw1.net
binomにある信頼区間をグラフ表示させてみた。
binom.ci <- function(x,n,cl=0.95){
PEci=binom::binom.confint(x,n,conf.level=cl)
y=nrow(PEci):1
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mar=c(3,3,2,2))
par(oma=c(0,4,0,0))
plot(PEci$upper,y,xlim=c(min(PEci$lower),max(PEci$upper)),type="n",
yaxt="n",ylab="",xlab="Confidence Interval",
main=paste("C.I. for", x,"out of",n))
points(PEci$upper,y,pch="|") ; points(PEci$lower,y,pch="|")
segments(PEci$lower,y,PEci$upper,y,lwd=3,
col=sample(colours(),1))
abline(v=PEci$mean[1],lty=3)
axis(2, at=y, labels=PEci$method,as.vector,las=2,cex=0.75)
PEci
}
binom.ci(3,312)
869:１３２人目の素数さん
20/05/08 00:48:23 CfFk/Uw1.net
>>826
情報、ありがとうございます。
?paletteのヘルプにあるスクリプトを実行すると体感できました。
#
870:# Demonstrate the colors 1:8 in different palettes using a custom matplot() sinplot <- function(main=NULL) { x <- outer( seq(-pi, pi, length.out = 50), seq(0, pi, length.out = 8), function(x, y) sin(x - y) ) matplot(x, type = "l", lwd = 4, lty = 1, col = 1:8, ylab = "", main=main) } sinplot("default palette") palette("R3"); sinplot("R3") palette("Okabe-Ito"); sinplot("Okabe-Ito") palette("Tableau") ; sinplot("Tableau") palette("default") # reset
871:１３２人目の素数さん
20/05/09 20:51:55 oDT7bFgO.net
rjagsのコードに
theta[i,j] <- min(1, exp(-alpha[i]*t[j]) + beta[i])といいれて、確率を１以下にして二項分布させてようとすると
Node inconsistent with parents というエラーがでた。 stackoverflow.comも明確な答が見いだせず、あれこれいじって
どうも確率1や0で二項分布する可能性があるとエラーがでるみたい、というのを発見。
1を0.99999とするとエラーが回避できた。
1を 1 - 1e-16としてもエラー回避できた。
872:１３２人目の素数さん
20/05/18 08:52:51.27 s/hrJxRo.net
mathematicaは内部どうなってんだか知らないが第何種ベッセル関数とか入ってくるめちゃくちゃ重い積分を計算してくれた
873:１３２人目の素数さん
20/05/25 12:44:17 3k8+9x6G.net
"
ある人物Dが新型コロナ肺炎に罹患したとする。
行動調査によって発症前にキャバクラに行っており
接客したキャバ嬢Cが人物D発症の２日後に発症していたことがわかった。
キャバ嬢Cは人物Dから移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢Cから移された可能性もあると主張して
その確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
"
潜伏期間が従う分布はこれ
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612
Gt <- function(delay){
C=rlnorm(1e6,ln_par1,ln_par2)
D=rlnorm(1e6,ln_par1,ln_par2)
mean(C-D > delay)
}
Gt(2)
874:１３２人目の素数さん
20/05/25 12:48:45 3k8+9x6G.net
怖い本で紹介されていたパッケージ polspline は秀逸。
発生させた乱数から確率密度関数を作ってくれる。
複雑な分布だと収束できず作れないこともあるけど
対数正規分布の差の分布を乱数発生させてpdfが作れた。
library(polspline)
delay=2
ln_par1 = 1.434065
ln_par2 = 0.6612
c=rlnorm(1e5,ln_par1,ln_par2)
d=rlnorm(1e5,ln_par1,ln_par2)
hist(c-d, freq=F,col='white',breaks=100,ylim=c(0,0.11),
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
fit=logspline(c-d)
curve(dlogspline(x,fit),-30,30,ann=F,bty='n',add=T)
segments(delay,0,delay,dlogspline(delay,fit),pch=19,col=4)
curve(dlogspline(x,fit),delay,30,add=T,type='h',col=4)
curve(plogspline(x,fit),-30,30,lty=3,bty='n')
1-plogspline(delay,fit)
fn <- function(delay) 1- plogspline(delay,fit)
curve(fn(x),0,14, bty='n' ,xlab='Delay', ylab='Probability')
875:１３２人目の素数さん
20/06/01 16:49:09 smFiMdov.net
今気づいたけどdplyr 1.0.0がcranに来てた
876:１３２人目の素数さん
20/07/25 08:42:04 1FL39oBq.net
Rのグラフの出力をGIF動画にしたかったので検索したら、このページが役立った。
URLﾘﾝｸ(cse.naro.affrc.go.jp)
877:ezawa/r-tips/r/36.html mytest <- function() { ANS <- readline("1+2? ") if (ANS == "3") cat("Correct!\n") else cat("Wrong...\n") } mytest() mymessage <- function() { Sys.sleep(3) cat("Hello.\n") } mymessage() 尚、GIF作製には https://www.screentogif.com/ がお手軽だった。
878:１３２人目の素数さん
20/07/26 10:34:51 FqjSGBlb.net
他人のコードを読んでいたら、
'%&%' <- function(x,y) paste0(x,y)
という関数が定義されていた。
これは便利と思った。
> '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi
[1] "√2 =1.4142135623731, π = 3.14159265358979"
879:１３２人目の素数さん
20/07/27 13:23:47 3h6isdRg.net
文字列をパイプでつないでいくのは面白いかも
880:１３２人目の素数さん
20/07/28 04:21:24 jhN+xo4C.net
'%.%' <- function(x,y) sum(x,y)
でベクトルの内積とかも
881:１３２人目の素数さん
20/07/28 04:22:11 jhN+xo4C.net
訂正
'%.%' <- function(x,y) sum(x*y)
882:１３２人目の素数さん
20/07/29 12:22:51 3HY+hEHZ.net
面白いけど使いどころが難しそう
%>%パイプとの親和性が無いのが痛い
883:１３２人目の素数さん
20/07/29 13:52:08 5tjjlndH.net
%文字%は二項演算子だから、定義次第じゃない？
884:１３２人目の素数さん
20/07/29 18:45:33 3HY+hEHZ.net
tidyverseや%>%を使わないならまあ
885:１３２人目の素数さん
20/07/30 18:02:20.97 8iPZt0Nw.net
>>835
それって意味があるの？
> paste0('√2 =', sqrt(2), ', π = ', pi)
[1] "√2 =1.4142135623731, π = 3.14159265358979"
普通にpaste0のまま使った方がタイプ量が圧倒的に少ないが。
886:１３２人目の素数さん
20/07/30 20:10:10 417El1m4.net
>>842
括弧の対応で混乱しそうなときにちょっと役にたつかも。
こんな感じ
eval(str2lang("expand.grid(" %&% paste0(rep("1:2",6),collapse=',') %&% ")"))
887:１３２人目の素数さん
20/07/31 07:51:08 0a0zx1wG.net
>>843
うーん、その例だとイマイチに思える
rerun(6, 1:2) %>% expand.grid
で済むものをわざわざ複雑にしているだけに感じる
888:１３２人目の素数さん
20/07/31 13:51:36.15 ImynyFWQ.net
>>844
無理やり、例示した感は否めないな。
標準装備でなら、
expand.grid(replicate(6,1:2,simplify = FALSE))
ですむし。
889:１３２人目の素数さん
20/07/31 23:36:09.03 0a0zx1wG.net
paste()やstr_c()以外を使うなら、sprintf()かglue::glue()が候補か
特にカンマが文字列に含まれるときにはpaste()より見やすくなるはず
sprintf("√2 = %s, π = %s", sqrt(2), pi)
glue("√2 = {sqrt(2)}, π = {pi}")
890:１３２人目の素数さん
20/08/01 18:51:48 2TFdzAyg.net
>>846
確かに便利ですね。
カンマを ','で入力するのは混乱のもとでしたから。
有用な情報、ありがとうございます。'
> glue::glue("√2 = {sqrt(2)}, sin(π/3) = {sin(pi/3)}")
√2 = 1.4142135623731, sin(π/3) = 0.866025403784439
891:１３２人目の素数さん
20/08/06 00:01:24 rivMs4GJ.net
初心者質問で申し訳ありません。
ある関数を作りたいのですが、その関数に入れる数値の数はランダムな場合はどのように書けばいいのでしょうか。
例えば、入れた要素をすべて掛け合わせる関数が作りたい場合、そこに入れる要素が（２，３，４）でも、（２，３，４，５，６）でも反応するような関数の
892:書き方を教えてくださると助かります。
893:１３２人目の素数さん
20/08/06 00:05:38 ID3ryEgK.net
ベクトルで渡せばいいのでは？
そういう話じゃない？
894:１３２人目の素数さん
20/08/06 01:18:04 rivMs4GJ.net
874さん。
847さん、回答ありがとうございます。
確かにベクトルで何かするのはわかるのですが、その部分がわからないのです。
出来れば下に問題を書くので、コードを書いていただけると非常に助かります。
問、正の整数であるベクトルを入力すると、そのベクトルの要素をすべて掛け合わせた値を結果として返す関数を作りなさい。
895:１３２人目の素数さん
20/08/06 02:03:47 V+bu3PeG.net
>>850
問題文これだけ？他に条件はないの？
標準の関数にあるけどいいのかな？
関数fを作るとして、一番いいのは標準の関数に余計な変更を加えないことだから、これが答えだけど…
f <- prod
896:１３２人目の素数さん
20/08/06 02:11:09 rivMs4GJ.net
849さん、返信ありがとうございます。
実は問２がありまして、たぶんそれは既存の関数にはないと思うので、そちらも教えていただけると非常に助かります。
問２、要素がすべて整数であるベクトルを入力すると、そのベクトルの要素の絶対値をすべて掛け合わせた値を結果として返す関数を作成せよ。
897:１３２人目の素数さん
20/08/06 05:13:12 03BiGX3I.net
>>852
f = function (x) abs(prod(x))
898:１３２人目の素数さん
20/08/06 05:15:08 03BiGX3I.net
>>851
while やfor もしくは再帰関数を使って作らせる問題かな？
899:１３２人目の素数さん
20/08/06 05:27:00.71 meNNWVIo.net
>>848
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans
}
> f(2,3,4)
[1] 24
> f(2,3,4,5,6)
[1] 720
900:１３２人目の素数さん
20/08/06 05:32:55.68 meNNWVIo.net
>>852
f <- function(...){
v=c(...)
ans=1
for(i in v) ans=ans*i
ans=ifelse(ans>0,ans,-ans)
ans
}
> f(-1,2,3)
[1] 6
> f(1 -2,-3,-4,-5)
[1] 60
901:１３２人目の素数さん
20/08/06 05:34:50.08 meNNWVIo.net
>>852
入力はベクトルだったから、
> f <- function(v){
+ ans=1
+ for(i in v) ans=ans*i
+ ans=ifelse(ans>0,ans,-ans)
+ ans
+ }
> f(c(-1,2,3))
[1] 6
> f(c(1 -2,-3,-4,-5))
[1] 60
902:１３２人目の素数さん
20/08/06 06:04:14.18 meNNWVIo.net
再帰だとこんな感じかな？
f <- function(...){
v=c(...)
n=length(v)
sub<-function(n){
if(n==0) return(1)
else v[n]*Recall(n-1)
}
sub(n)
}
> f(2,3,4)
[1] 24
> f(2,3,4,5,6)
[1] 720
903:１３２人目の素数さん
20/08/06 06:35:53 eFxPvNIl.net
みんなアタマいいなぁ…
904:１３２人目の素数さん
20/08/06 08:34:19 E/ufUouA.net
>>859
854, 855, 856を真似るなよ。
905:１３２人目の素数さん
20/08/06 08:49:09 E/ufUouA.net
再帰版はこんな感じかな。
rec.prod <- function(x) {
# 再帰により x の要素の絶対値の積を求める。
# 実際には prod(abs(x)) が良い。
l <- length(x) # ベクトルの要素数
if (l == 0) # 空のベクトルなら
return(1) # 積は1
x1.a <- abs(x[1])
if (l == 1)
x1.a
else
x1.a * rec.prod(x[-1])
}
906:１３２人目の素数さん
20/08/06 11:42:21 rivMs4GJ.net
質問者本人です。
皆様のお力添えのおかげで無事問題を解くことができました。
初心者質問にやさしく回答して頂いて感謝しかありません。
本当にありがとうございました。
907:１３２人目の素数さん
20/08/06 19:28:54.17 meNNWVIo.net
>>861
この方が>858みたいに関数の中で関数を定義とかなくてカッコいいと思ったのでsystem.timeで時間を比較しようと思って実験を始めたら、
何故か、rec.prodの方がoverflowしやすいなぁ。
頭からかけても尻からかけても同じじゃないかと思うんだが。
> prod(1:100)
[1] 9.332622e+157
> rec.prod <- function(x) {
+ # 再帰により x の要素の絶対値の積を求める。
+ # 実際には prod(abs(x)) が良い。
+ l <- length(x) # ベクトルの要素数
+ if (l == 0) # 空のベクトルなら
+ return(1) # 積は1
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rec.prod(x[-1])
+ }
> rec.prod(1:100)
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
> f <- function(...){
+ v=c(...)
+ n=length(v)
+ sub<-function(n){
+ if(n==0) return(1)
+ else v[n]*Recall(n-1)
+ }
+ sub(n)
+ }
> f(1:100)
[1] 9.332622e+157
908:１３２人目の素数さん
20/08/06 20:06:06 meNNWVIo.net
比較のために、関数名の�
909:ｷさを同じにして、時間のかかるRecallを関数名にしてsystem.timeで比較してみた。 > rm(list=ls()) > rp <- function(x) { + # 再帰により x の要素の絶対値の積を求める。 + # 実際には prod(abs(x)) が良い。 + l <- length(x) # ベクトルの要素数 + if (l == 0) # 空のベクトルなら + return(1) # 積は1 + x1.a <- abs(x[1]) + if (l == 1) + x1.a + else + x1.a * rp(x[-1]) + } > f1 <- function(...){ + v=c(...) + n=length(v) + sub<-function(n){ + if(n==0) return(1) + else v[n]*Recall(n-1) + } + sub(n) + } > f2 <- function(...){ + v=c(...) + n=length(v) + sub<-function(n){ + if(n==0) return(1) + else v[n]*sub(n-1) + } + sub(n) + } > system.time(replicate(1e5,prod(1:10))) user system elapsed 0.11 0.00 0.11 > system.time(replicate(1e5,rp(1:10))) user system elapsed 2.64 0.00 2.64 > system.time(replicate(1e5,f1(1:10))) user system elapsed 2.44 0.01 2.47 > system.time(replicate(1e5,f2(1:10))) user system elapsed 1.86 0.00 1.87
910:１３２人目の素数さん
20/08/06 20:31:17.99 E/ufUouA.net
それはね、856は厳密な意味での再帰ではないから速いしリ、ソースを食わないんだよ。ループを再帰風に書いただけだから。それよりもループのほうが速いに決まっている。
この問題を再帰で解くのはムダ以外の何ものでもない。
911:１３２人目の素数さん
20/08/06 20:33:25.51 w8yiRRv5.net
>>865
喩えれば、
親分が子分を再帰で酷使しているだけだから親分は消耗しない
という理解でいい？
912:１３２人目の素数さん
20/08/06 20:47:18.11 E/ufUouA.net
全然違う。
再帰は、複雑な問題をより小さい問題に分割することがキモだ。
856は元のベクトルはそのまま変わらずに存在していて、ある特定の要素を一つずつ処理しているからループとやってることは変わりない。
859は、一つ短いベクトルを自分自身に処理させているから真の再帰なのだ。
913:１３２人目の素数さん
20/08/06 22:25:29 meNNWVIo.net
>>867
f自身は再帰関数じゃないけど、呼び出しているsubは再帰関数では？
こういうのと同じ。
f <- function(x){ # 数列の和の階乗を返す
m=sum(x)
sub <- function(n){
if(n==1) return(1)
else{
return(n*Recall(n-1))
}
}
fact(m)
}
> f(c(1,2,3))
[1] 720
914:１３２人目の素数さん
20/08/06 22:50:49 XP8VEZbb.net
>>867
(ご入力訂正)
f自身は再帰関数じゃないけど、呼び出しているsubは再帰関数では？
こういうのと同じ。
f <- function(x){ # 数列の和の階乗を返す
m=sum(x)
sub <- function(n){
if(n==1) return(1)
else{
return(n*Recall(n-1))
}
}
sub(m)
}
> f(c(1,2,3))
[1] 720
915:１３２人目の素数さん
20/08/07 21:08:59 S7qsZ31k.net
856はsub()の引数がvでなくnなことと、nの入力に対して関数外のv[n]が帰ってくるのが、少し気になるかも
856が再帰っぽくないってのは、859はxを分割・再帰していってるのに対して、856はnを減じていってその都度vのn番目の値を参照してるところかな
それから、859がrec.prod(1:13)であっさりオーバーフローしてしまうのは、integer型だから
あと実行速度が少し遅いのは、xl.aへの代入操作と、if(l ==1){}の余分な条件分岐のせいだね
その辺を修正してやれば、こんな感じかな
# 再帰で問２
f <- function(x){
if(length(x) == 0) return(1)
abs(x[1]) * f(x[-1])
}
916:１３２人目の素数さん
20/08/08 05:40:19.06 2ggSSq05.net
>>870
勉強になりました。
ついでにwhile版を書いてみました。
f <- function(x){
ans=1
while(length(x)){
ans=ans*abs(x[1])
x=x[-1]
}
ans
}
917:１３２人目の素数さん
20/08/08 09:18:18 FvlAeRBY.net
>>870
859が遅い根本の原因は、x[-1] が生成されること。
他の指摘は、細かい高速化の工夫だから初心者に薦めるのはどうかと。
現に869で全く実用的に意味のないコードを書いているし。
題意を素直に表したほうがよいのでは？
絶対値の積は、積の絶対値と等しいから、absを呼ぶのは1回で済むが、何も説明無しにそうするのは私的にはアウトだ。メンテナンスの時に自分や他人が苦しむから。
918:１３２人目の素数さん
20/08/08 11:01:30 559rm1Tc.net
>>872
意味不明なこと言いだしてるが大丈夫か、落ち着け
859は、無駄な条件分岐あるし、すぐオーバーフローするし、コードが冗長だしで、それこそ初心者以外にもNGだぞ
絶対値の積うんぬんに関しては、859と同じでabsを毎回呼んでる、余計な代入操作が無いだけだ
落ち着いてコードを見ろ
速度うんぬんの話はこれまで話題に出てきたように、856と比べての話だ、話についてこれないなら無理に入ってこなくていい
別に高速化を目的としてるわけじゃなく、859から冗長な部分やオーバーフローの原因を取り除いた結果高速化されたってだけだ
だが、速度が遅い原因がxl.aへの代入でなく、x[-1]作成が原因てのはそのとおりだ
文章を書いてる途中で、オーバーフローの原因（代入操作はこっちだ）と勘違いしてしまったようだ
これは訂正に感謝だ
919:１３２人目の素数さん
20/08/08 11:24:32 FvlAeRBY.net
>>873
おまえこそもちつけ。
再帰は、再帰の停止条件と、自分自身が行う作業に分離できる。
この場合の停止条件は、要素が1つになったときと考えるのが自然だろう。
数学的には要素が0でもよいし、Rの場合はそれでうまくいくから問題はない。だが、そういった詳細を何も説明せずにこれでよい、と示すのはどうかな？
あと、絶対値を取ったのを代入しているのは、この作業が複雑になった場合の修正の確実さを優先するためで、優秀なバイトコンパイラであれば速度の差はなくなる。
920:１３２人目の素数さん
20/08/08 23:06:27.95 UlZ+DmHh.net
一人だけエラー吐いて、まともに動かない関数書いた無能が、
それでも、恥ずかしげもなく、一生懸命マウント取ろうとあがく姿は、
滑稽で哀れになってくるわ
921:１３２人目の素数さん
20/08/13 17:34:36 5b54QocS.net
＞859がrec.prod(1:13)であっさりオーバーフローしてしまうのは、integer型だから
の意味が理解できなくて、実験してみた。
> rec.prod <- function(x) {
+ l <- length(x)
+ if (l == 0)
+ return(1)
+ x1.a <- abs(x[1])
+ if (l == 1)
+ x1.a
+ else
+ x1.a * rec.prod(x[-1])
+ }
>
>
> rec.prod(1:13)
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
> rec.prod(as.numeric(1:13))
[1] 6227020800
> rec.prod(c(1,2,3,4,5,6,7,8,9,10,11,12,13))
[1] 6227020800
> rec.prod(c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L))
[1] NA
Warning message:
In x1.a * rec.prod(x[-1]) : NAs produced by integer overflow
integer型のデータだとオーバーフローしちゃうけど、as.numericでdouble型にするとエラーがでない。
一方、
> f <- function(x){
+ if(length(x) == 0) return(1)
+ abs(x[1]) * f(x[-1])
+ }
> f(1:13)
[1] 6227020800
> f(as.numeric(1:13))
[1] 6227020800
> f(c(1,2,3,4,5,6,7,8,9,10,11,12,13))
[1] 6227020800
> f(c(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L))
[1] 6227020800
>
とエラーがでないなぁ。
何故だろう？
922:１３２人目の素数さん
20/08/20 13:26:56 mPUCNi08.net
Package ‘XXX’ was installed before R 4.0.0: please re-install it
とメッセージがでてXXXをインストールすると次は
Package ‘YYY’ was installed before R 4.0.0: please re-install it
とでてきてうんざりしていた。
stackoverflow.comで以下の答をみつけて解決した。
.libPaths()
update.packages(ask=FALSE, checkBuilt=TRUE)
923:１３２人目の素数さん
20/09/09 22:54:07 IR7822fG.net
ありがとうございます
924:１３２人目の素数さん
20/10/19 07:15:12.45 21zikI9w.net
PならばQ をパイプで P %=>% Q として論理演算で遊んでみた。
"
問題 ： 「シリツ医 ならば (馬鹿 ならば 裏口 である)」という命題と同値な命題はどれか？
1 : シリツ医 ならば (裏口 ならば 馬鹿 である)
2 : 馬鹿 ならば (シリツ医 ならば 裏口 である)
3 : 馬鹿 ならば (裏口 ならば シリツ医 である)
4 : 裏口 ならば (シリツ医 ならば 馬鹿 である)
5 : 裏口 ならば (馬鹿 ならば シリツ医 である)
"
# PならばQ ≡ (P かつ (Qでない))ではない
'%=>%' = function(P,Q) !(P & !Q)
# premise
A = function(S,
925:B,U) S %=>% (B %=>% U) # choice B1 = function(S,B,U) S %=>% (U %=>% B) B2 = function(S,B,U) B %=>% (S %=>% U) B3 = function(S,B,U) B %=>% (U %=>% S) B4 = function(S,B,U) U %=>% (S %=>% B) B5 = function(S,B,U) U %=>% (B %=>% S) # premise => choice C1 = function(S,B,U) A(S,B,U) %=>% B1(S,B,U) C2 = function(S,B,U) A(S,B,U) %=>% B2(S,B,U) C3 = function(S,B,U) A(S,B,U) %=>% B3(S,B,U) C4 = function(S,B,U) A(S,B,U) %=>% B4 (S,B,U) C5 = function(S,B,U) A(S,B,U) %=>% B5 (S,B,U) # combination grid of TRUE & FALSE gr=expand.grid(c(T,F),c(T,F),c(T,F)) # test for premise => choice all(mapply(C1,gr[,1],gr[,2],gr[,3])) all(mapply(C2,gr[,1],gr[,2],gr[,3])) all(mapply(C3,gr[,1],gr[,2],gr[,3])) all(mapply(C4,gr[,1],gr[,2],gr[,3])) all(mapply(C5,gr[,1],gr[,2],gr[,3]))
926:１３２人目の素数さん
20/10/19 07:15:28.30 21zikI9w.net
# choice -> premise
D1 = function(S,B,U) B1(S,B,U) %=>% A(S,B,U)
D2 = function(S,B,U) B2(S,B,U) %=>% A(S,B,U)
D3 = function(S,B,U) B3(S,B,U) %=>% A(S,B,U)
D4 = function(S,B,U) B4(S,B,U) %=>% A(S,B,U)
D5 = function(S,B,U) B5(S,B,U) %=>% A(S,B,U)
# test for choice -> premise
all(mapply(D1,gr[,1],gr[,2],gr[,3]))
all(mapply(D2,gr[,1],gr[,2],gr[,3]))
all(mapply(D3,gr[,1],gr[,2],gr[,3]))
all(mapply(D4,gr[,1],gr[,2],gr[,3]))
all(mapply(D5,gr[,1],gr[,2],gr[,3]))
927:１３２人目の素数さん
20/11/17 19:55:57.27 UJMPy762.net
,で文字列を分けようとしたらうまくいかなかったがHELPをみたら解説してあった。
> #.で文字列を分ける
> unlist(strsplit('3.14','.'))
[1] "" "" "" ""
> unlist(strsplit('3.14','\.')
Error: '\.' is an unrecognized escape in character string starting "'\."
>
> unlist(strsplit('3.14','[.]'))
[1] "3" "14"
> unlist(strsplit('3.14','.',fixed=TRUE))
[1] "3" "14"
928:１３２人目の素数さん
20/11/17 20:18:32.66 /TSutyYc.net
unlist(strsplit("3.14", "¥¥."))
人に見てもらいたかったら読みやすく書けよ。アンタのプログラムはいつも読みにくい。
929:１３２人目の素数さん
20/11/18 01:42:18.74 CJUGi+QQ.net
>>882
備忘録として書いたのでレスがつくとは思ってなかったが、レスthanx.
930:１３２人目の素数さん
20/12/10 02:17:40.15 Whv7aYu3.net
数値積分をさせていたら
初めてみるこんなエラーが返ってきた。
roundoff error is detected in the extrapolation table
困ったときのstackoverflow.comで
URLﾘﾝｸ(stackoverflow.com)
As for the integration problem, package cubature generally does a better job when integrate
と解決法が記載されていた
integrate(fn, lower, upper)$value
でエラーが返ってきたが、
cubintegrate(fn, lower,upper, method='pcubature')$integral
とするとエラーなしに計算してくれた。
（備忘録）
931:１３２人目の素数さん
20/12/12 23:25:15.53 DvLXeu4z.net
それエラーが出るべきなのに出ない誤りなんじゃ、って可能性はないの？
932:１３２人目の素数さん
20/12/13 08:12:33.91 xmsEH+j9.net
>>885
別スレの問題を「プログラムで面積を数値計算させていたらエラー発生
体験したければ読みにいコードの下記をどうぞ。
r= 0.66017210648907 ; s=2.89088405538017
# 面積計算にintegrate で エラー
sub <- function(c,Area1=1){ # c 円Cの中心のｘ座標
b0=(c^2+r^2-s^2)/(2*c) # B(b0,b1) # 円Cと円Oの交点Bの座標
b1=sqrt(r^2-(c^2+r^2-s^2)^2/(4*c^2))
blue = function(x) sqrt(s^2-(x-c)^2) # 円Cの方程式
BLUE= integrate(blue,b0,c+s)$value*2 # 青の面積
green = function(x) sqrt(r^2-x^2) # 円Oの方程式
GRE
933:EN=integrate(green,-r,b0)$value*2 # 緑の面積 GREEN+BLUE-Area1 # GREEN+BLUE=1 } sub=Vectorize(sub) x=seq(-r-s,r-s,length=100) ; plot(x,sub(x),pch=19) ; abline(h=0,lty=3) # 面積計算にcubintegrate でエラーなし library(cubature) sub <- function(c,Area1=1){ # c 円Cの中心のｘ座標 b0=(c^2+r^2-s^2)/(2*c) # B(b0,b1) # 円Cと円Oの交点Bの座標 b1=sqrt(r^2-(c^2+r^2-s^2)^2/(4*c^2)) blue = function(x) sqrt(s^2-(x-c)^2) # 円Cの方程式 BLUE= cubintegrate(blue,b0,c+s,method="pcubature")$integral*2 # 青の面積 green = function(x) sqrt(r^2-x^2) # 円Oの方程式 GREEN=cubintegrate(green,-r,b0,method='pcubature')$integral*2 # 緑の面積 GREEN+BLUE-Area1 # GREEN+BLUE=1 } sub=Vectorize(sub) x=seq(-r-s,r-s,length=100) ; plot(x,sub(x),pch=19) ; abline(h=0,lty=3)
934:１３２人目の素数さん
20/12/13 08:15:08.06 xmsEH+j9.net
入力の省力化にこんな関数を道具箱に詰めた
# 定積分
integral <- function(fn,lwr,upr,...){
cubature::cubintegrate(fn,lwr,upr,method='pcubature',...)$integral
}
935:１３２人目の素数さん
20/12/13 09:32:37.11 xmsEH+j9.net
>>886
library(MASS)にareaという面積計算の関数があるのを思い出したのでやってみた。
library(MASS)
r= 0.66017210648907 ; s=2.89088405538017
# 面積計算にMASS::area
sub <- function(c,Area1=1){ # c 円Cの中心のｘ座標
b0=(c^2+r^2-s^2)/(2*c) # B(b0,b1) # 円Cと円Oの交点Bの座標
b1=sqrt(r^2-(c^2+r^2-s^2)^2/(4*c^2))
blue = function(x) sqrt(s^2-(x-c)^2) # 円Cの方程式
BLUE=area(blue,b0,c+s)*2 # 青の面積
green = function(x) sqrt(r^2-x^2) # 円Oの方程式
GREEN=area(green,-r,b0)*2 # 緑の面積
GREEN+BLUE-Area1 # GREEN+BLUE=1
}
sub=Vectorize(sub)
x=seq(-r-s,r-s,length=100) ; plot(x,sub(x),pch=19) ; abline(h=0,lty=3)
こっちはエラーがでない。
936:１３２人目の素数さん
20/12/21 20:43:26.34 gcvC1Sv7.net
func <- function (...) {
return(anova(...))
}
func(fit1, fit2, fit3)
オブジェクトの引数をうけとって、関数内の関数anovaに渡して
ちゃんと動かすにはどうしたらいいですか？
937:１３２人目の素数さん
20/12/23 07:21:19.77 N7Fbsys3.net
>>889
>775みたいに
c(...)もしくはlist(...)で引き出すのでは動作しない？
938:１３２人目の素数さん
20/12/23 11:45:33.99 DPeLDm1M.net
>>889
> mydata <- data.frame(x = runif(10), y = runif(10), z = runif(10))
> xy <- lm(y ~ x, data = mydata)
> zy <- lm(y ~ z, data = mydata)
> func <- function (...) {
return(anova(...))
}
> func(xy, zy)
Analysis of Variance Table
Model 1: y ~ x
Model 2: y ~ z
Res.Df RSS Df Sum of Sq F Pr(>F)
1 8 0.40535
2 8 0.39312 0 0.012234
ちゃんと動くけど。
何がどう問題なのかを明確にして質問しないと助言を得られないよ。
939:１３２人目の素数さん
21/01/03 22:00:14.73 8tLYm46h.net
やっぱり、パッケージになっている関数は高速だな。
1万までの素数を求める
> rm(list=ls())
> primes <- function(n) (1:n)[-outer(2:n,2:n)][-1]
> system.time(primes(10000))
user system elapsed
0.60 1.92 2.88
> rm(list=ls())
> library(numbers)
> system.time(Primes(10000))
user system elapsed
0.00 0.00 0.03
940:１３２人目の素数さん
21/01/21 00:00:54.88 KZXQ/VTY.net
時刻と、その時に従業員がいた場所（住所）がデータとして一覧化されていて、
いつ、どこによくいるかをクラスタリングしたいのですがRは出来るのでしょうか？
941:１３２人目の素数さん
21/01/24 12:08:49.15 oqt0Qph1.net
>>893
SpaTimeClusパッケージ
942:１３２人目の素数さん
21/02/14 11:56:55.00 g9kE6hlo.net
IPB size
1 I 82571
2 B 4909
3 B 230
4 B 231
こんな感じのデータフレームを複数ファイルを読み込んで
ファイル別に横につなげる場合どのようにfor文を書けばいいですか?
943:１３２人目の素数さん
21/02/14 12:28:29.31 n5qxzAHa.net
purrrパッケージを使った方が速いんじゃね？
944:１３２人目の素数さん
21/02/14 16:02:48.85 g9kE6hlo.net
> 894
ありがとう。これで無事横並べられた
URLﾘﾝｸ(gist.github.com)
945:１３２人目の素数さん
21/03/16 13:05:28.85 cJpJJpsg.net
# 振幅１の正弦波の一周期の長さを求めよ
の計算に、一般化して計算できるように関数を書いてみたのだが、
文字列での入力が必要でどうも美しくない。
# 関数y=f(x)の曲線でx=from からx=toまでの長さを求める
CurveLength <- function(f='sin(x)',from=-pi,to=pi){
str=paste("deriv(~ ",f,",","'x',func=TRUE)")
Df <- eval(str2lang(str))
f1 <- function(x) as.numeric(attributes(Df(x)))
f1=Vectorize(f1)
integrate(function(x) sqrt(1+f1(x)^2), from, to, rel.tol = 1e-12)
}
計算としてはあっていると思う。
> CurveLength()
7.640396 with absolute error < 2.6e-13
> CurveLength('log(1+cos(x))',0,pi/2)
1.762747 with absolute error < 2e-14
もっとエレガントなコードがあればご教示をお願いします。
946:１３２人目の素数さん
21/03/18 08:46:13.56 IPPbuENC.net
時系列データ成型の解説でいいサイトないですか
947:１３２人目の素数さん
21/03/31 01:41:06.70 Cwg7NuDW.net
>>898
コードのことならプログラム板のほうで訊いたほうがよさそう
948:１３２人目の素数さん
21/06/05 21:13:05.02 5b0TiKZo.net
R version 4.1.0 (2021-05-18) -- "Camp Pontanezen" になってパイプに標準で対応したのは嬉しい。
function(x) が \(x)と略記可能になった。
ちょっと練習
f <- \(x) sample(100,x,rep=TRUE)
f(50) |> sort() |> unique() |> length()
> f <- \(x) sample(100,x,rep=TRUE)
> f(50) |> sort() |> unique() |> length()
[1] 40
詳細は
URLﾘﾝｸ(cran.r-project.org)
で
949:１３２人目の素数さん
21/06/06 12:41:44.43 X+DUyz8T.net
Vectorize()もパイプで受けてくれた。
(\(n) sample(10,n, replace = TRUE) |> unique() |> sort()) |> Vectorize() -> f.new
f.new(1:5)
f.old <- Vectorize(function(n){
sort(unique(sample(10,n,replace = TRUE)))})
f.old(1:5)
950:１３２人目の素数さん
21/06/06 12:54:15.47 X+DUyz8T.net
|> curve()は描画されなかった。
(\(x) sqrt(1-x^2)) |> Vectorize() |> curve()
> (\(x) sqrt(1-x^2)) |> Vectorize()) |> curve()
Error: unexpected ')' in "(\(x) sqrt(1-x^2)) |> Vectorize())"
(\(x) sqrt(1-x^2)) |> Vectorize() -> fn
curve(fn)
1/4円を描画
951:１３２人目の素数さん
21/06/06 22:14:15.22 X+DUyz8T.net
パイプで速度が落ちるかと思ったけど変わらないね。
> (\(x) sample(1000,x,replace=TRUE) |> sort() |> unique() |> length()) -> f.new
> system.time(replicate(1e5,f.new(5000)))
user system elapsed
28.66 0.01 28.70
> f.old <- function(x) length(unique(sort(sample(1000,x,replace=TRUE))))
> system.time(replicate(1e5,f.old(5000)))
user system elapsed
28.65 0.00 28.67
952:１３２人目の素数さん
21/06/06 22:37:59.70 jjPG368x.net
そのパイプを使ったコードをquote関数に噛ませば分かるけど実際はパーサーが下のコードに変換してるだけらしい。
なので実態としては演算子や関数ではないのでオーバーヘッドが生じないとのこと。
953:１３２人目の素数さん
21/06/08 05:41:54.05 uSV089kB.net
>>905
レスありがとうございます。早速、やってみました。
> quote((\(n) sample(10,n, replace = TRUE) |> unique() |> sort()) |> Vectorize())
Vectorize((function(n) sort(unique(sample(10, n, replace = TRUE)))))
954:１３２人目の素数さん
21/06/08 06:02:44.83 uSV089kB.net
dplyrのパイプ%>%だと他に引数がないときは関数の()は省略できるけど
今回、搭載された|>は()がないとエラーになるな。
> library(dplyr)
> sample(10,10,replace=TRUE) %>% sort
[1] 1 2 4 5 6 7 7 7 7 8
> sample(10,10,replace=TRUE) %>% sort(decreasing = TRUE)
[1] 9 9 7 6 5 5 3 1 1 1
> sample(10,10,replace=TRUE) |> sort
Error: The pipe operator requires a function call as RHS
> sample(10,10,replace=TRUE) |> sort()
[1] 3 4 5 5 6 6 8 8 8 10
> sample(10,10,replace=TRUE) |> sort(decreasing = TRUE)
[1] 10 10 9 8 7 7 6 6 5 2
955:１３２人目の素数さん
21/06/08 06:10:22.55 uSV089kB.net
パイプで送る関数に引数が１個なら自分でパイプを定義するのと変わらないなぁ。
オンラインで実行できるサイトだと必ずしも最新版のRに対応していないし、外部ライブラリ非対応のところもあるので
自分でパイプを定義の方がいいかも
'%|%' <- function(x,FUN) FUN(x)
> sample(10,10,replace=TRUE) %|% sort %|% unique
[1] 1 3 5 6 7
956:１３２人目の素数さん
21/08/29 23:23:03.95 OFN46wh5.net
# 備忘録
hogehoge <- function(x){
cat(deparse(substitute(x)),'=',x,'\n')
}
ch='Swiss'
hogehoge(ch)
> hogehoge(ch)
ch = Swiss
957:１３２人目の素数さん
21/11/05 22:47:20.72 lTyg7nhh.net
これ↓、思ったように動かないけど、普通？
iris[any(iris[,3]==c(6.1,5.1)) ,]
いまいちR言語の仕様がわからん。
958:１３２人目の素数さん
21/11/05 23:38:14.16 I43ryX9p.net
何をしたいのか説明できない理系脳かな？
適当に推測するに、
iris[iris[, 3] %in% c(6.1, 5.1), ]
959:１３２人目の素数さん
21/11/06 21:41:29.28 8hjJFETA.net
ありがとう。
雰囲気で伝わるかなと思った。
960:１３２人目の素数さん
21/11/06 21:45:56.76 8hjJFETA.net
SQL(例えばAccess)とRの、データ分析に関する違いって何だと思う？
Rの有利な点は表を容易に分割→編集→結合できる点で、Accessの有利な点はクエリを何段にも重ねることで表を視覚的に扱えることだと思うんだけど、どうかな？
要約すると、細かいとこに手が届くのと、そうではないけど直感的なところ。
961:１３２人目の素数さん
21/11/06 21:53:39.63 8hjJFETA.net
SQLがAccessの例に限定されちゃってるね。
でも同じように文字で書くなら、Rの方が便利に思えるな。
Officeが使える環境でR使ってるから、Rだけの環境に移行すればどんな不便が発生するか検討中なんだけど、どう思う？
962:１３２人目の素数さん
21/11/08 19:14:28.35 xqrS4yV9.net
>>913
> SQL(例えばAccess)とRの、データ分析に関する違いって何だと思う？
何を言っているのか分からないが、用意されいている分析手法の豊富さが全く違うだろ。
例えば、多重検定の検定法とかSQLにあるの？ノンパラメトリックの多重検定とかちゃんとできるの？
963:１３２人目の素数さん
21/11/08 19:14:46.14 xqrS4yV9.net
>>913
> SQL(例えばAccess)とRの、データ分析に関する違いって何だと思う？
何を言っているのか分からないが、用意されいている分析手法の豊富さが全く違うだろ。
例えば、多重検定の検定法とかSQLにあるの？ノンパラメトリックの多重検定とかちゃんとできるの？
964:１３２人目の素数さん
21/11/09 20:40:06.49 hGmtIAbe.net
そうだよね。
組み合わせて使うもんだよね。SQLとRは。
ただ、境界線をどこで引くべきか微妙なところで悩む。RだとRStudioのレポートとして残せるから、分析手順を残していくことを考えると極力R使った方がいいのかな。
965:１３２人目の素数さん
21/12/04 13:16:48.05 3GLpEIlX.net
いや、SQLとRは全くの別物だろ
SQLはデータベースのマネジメントシステムで、
Rはデータ分析のツールじゃん
966:１３２人目の素数さん
21/12/04 14:11:32.93 0kKQqVbM.net
別物だから組み合わせて使う意味がある
967:１３２人目の素数さん
21/12/04 20:41:11.61 fCvhwSlW.net
いや、SQLとRは全くの別物だろ
組み合わせて使う意味はない
968:１３２人目の素数さん
21/12/04 22:13:49.92 6B5EP1qJ.net
そうだな、お前は使わなくていい
絶対に使うな
969:１３２人目の素数さん
21/12/04 22:16:19.95 sikGs8rn.net
別物だから別々に使う意味がある
970:１３２人目の素数さん
22/02/28 07:22:18.49 A7JZ2ass.net
URLﾘﾝｸ(twitter.com)
を実行出来た方いますか？
(deleted an unsolicited ad)
971:sage
22/02/28 12:38:08.27 dt49D0nX.net
URLﾘﾝｸ(detail.chiebukuro.yahoo.co.jp)
で解決しました。
スレを汚してすみませんせした。
972:１３２人目の素数さん
22/03/02 03:27:08.75 RV9q7yNr.net
今更レスだけどなんでaccessと比べるのか、ex
973:celじゃねーの excelは商用だけあって機能的には検定回帰成分分析因子分析何でもできてRに劣るところは無いし、統計以外なら最適化が強い Rより遅いしVBAは醜いし、365なんていきなり環境ぶっ壊されかねないからRに安住するけど
974:１３２人目の素数さん
22/03/02 03:30:48.43 RV9q7yNr.net
インタラクティブに使う分にはCSE使えるexcelの方が便利だな
ぜひともRに取り入れたい機能の一つ
975:１３２人目の素数さん
22/03/02 03:42:31.09 RV9q7yNr.net
excelで巨大データ扱うにはaccessインテグレーション機能で補完し合うといいとは聞くけど、俺のアカデミック割のofficeではaccessはハブられてた悲しみ
976:１３２人目の素数さん
22/03/06 16:53:40.50 BGFbnw+F.net
精度が怪しいだけ
977:１３２人目の素数さん
22/03/06 16:55:06.71 BGFbnw+F.net
コマンドで操作しづらい
978:１３２人目の素数さん
22/03/18 15:05:36.24 Z+B2NJF6.net
contourのzlimはどう設定すればよいのでしょうか？
URLﾘﾝｸ(navaclass.com)で
ユークリッド距離のzlimは
contour(x,y,z.E,zlim=c(0,3),nlevels = 5,add=T)
のようにzlim=c(0,3)に設定されているのですが、
マハラノビス距離のzlimは
contour(x,y,z.D,zlim=c(0,1),nlevels = 5,labels="", add=T)
のようにzlim=c(0,1)に設定されています。
確かに、ユークリッド距離のほうは(1, 1)の点が1.0の等高線上にありますし、
マハラノビス距離のほうもAは0.75、Bは2.2という距離に相応しいところに点があります。
ただ、zlim=c(0,3)の意味はきっと0から3までの範囲に限るという意味なんでしょうが、
z.Eとの値も見ても3で区切られているようには見えません。
実験してみましたが、z.Eのほうのzlimは3.5、z.Dのほうのzlimは1.4で等高線が変化するようです。
これでは新たな相関係数を使ったときにどう設定すればよいのか分かりません。
どうか教えて下さい。
979:１３２人目の素数さん
22/03/18 15:25:12.51 3f0eXSLb.net
>>930
xlim, ylimは指定した値でスパッと切るわけではないので（前後に余裕を持たせるようになっている）、
zlimも同じ仕様だと想像するよ。
xlimやylimで指定した特定の値でスパッと切るオプションがあったと思うので、
zlimにも有効でないかどうか、そこから調べて見れば？
980:１３２人目の素数さん
22/03/18 17:48:09.58 Z+B2NJF6.net
>>931
ありがとうございます。
今、いろいろ計算した結果、等高線の位置自体は正しいことが判りました。
ただ、問題は等高線の刻みがzlimと相関係数の組み合わせによって変わることです。
具体的には、等高線が0.5刻みになったり1.0刻みになったり0.25刻みになったりします。
にもかかわらずnlevelsは5本のままなので、全体が膨らんだり縮んだり、見た目が全然違ってくるんですよね。
その度に手でzlimの値を1とか3とか5に変えてあげないといけません…これらの数字はまったく理解できません。
なので質問を変えさせていただいて、
相関係数が0以外の時でも等高線上に0.5, 1.0, 1.5, ...などのラベルを強制的に表示する方法は無いでしょうか？
相関係数が0（つまりユークリッド距離）の時は等高線上にそれらの数値が表示されます。
それが相関係数が0.1にでもなった途端に表示されなくなります
（数値表示用に輪が途切れた形でブランクは空いているようなのですが）。
等高線上に0.5, 1.0, 1.5, ...などのラベルが表示されるようになれば、上記の問題も許容できると思います。
余談ですが、
>スパッと切るオプション
調べてみると、xaxs="i"のことのようですね。
"r"だとグラフの端に4%の余白があったのが"i"だと無くなりました。
981:１３２人目の素数さん
22/03/18 18:10:24.04 Z+B2NJF6.net
すみません、labels=“”になっているというオチでした。
とりあえず、これで問題ないか確認します。
982:１３２人目の素数さん
22/03/21 13:18:24.75 k7XqFrKG.net
ご助力をいただきたく、質問いたします。
対象データ
データフレームA内のB列, <chr>型 中身は例として 文字列で"2022/02/24"が入っている。
やりたいこと
この文字列をRのDate型に変換したい
やってみたこと。
tidyverseを用いて
A %>% select(B) %>% as.Date()
これだとエラーが出ました。
エラー文は「'.'からクラス"Date"へ変換はできません」
どうやったらこのB列をRのDate型に変換できますでしょうか。
983:１３２人目の素数さん
22/03/21 15:20:29.09 SKGsv+wa.net
>>934
パイプ処理したいならlubridateパッケージが便利です。チートシート（英語）が公開されているので、使い方は確認してみてください。
パッケージ名のスペル間違ってたらごめん。
984:１３２人目の素数さん
22/03/21 17:46:23.39 qn67+d1x.net
>>934
変換できない理由は、dplyr::select()関数の返り値がdata.frame型で
as.Date()関数の引数がベクトル型でなければならないからです。
なので、
as.Date(A$B)
とすれば、Date型でベクトルな返り値を得ることができます。tidyverseを
使いたい（というかパイプで処理したい）なら
A %>%
.$B %>% # ベクトルとして取り出す
as.Date()
となります。
返り値をdata.frame型で得たい場合は、dplyr::mutate()関数を使って
A %>%
dplyr::select(B) %>%
dplyr::mutate(B = as.Date(B))
とします。
tidyverseは、基本的にデータフレームを入力としデータフレームを出力と
しますので、列をベクトルとして処理したい場合はmutate()関数を使う必要が
あります。
985:１３２人目の素数さん
22/03/21 20:28:10.47 k7XqFrKG.net
>935様
>936様
ご教示ありがとうございます。
参考にまずは試してみます。
986:１３２人目の素数さん
22/03/22 07:49:54.09 kcYWMWlK.net
>>934
selectじゃなくてpullにすれば動くよ
987:１３２人目の素数さん
22/03/22 12:12:27.45 CNKDZNvw.net
pull()は知らなかったけど引数を二つ指定するとnamed vectorも作れるんだな。
988:１３２人目の素数さん
22/03/25 19:24:11.66 YaAd0nW/.net
借金を返す手段は大増か
踏み倒し（＝ハイパーインフレ）の２者選択。
今、日本はお金をどんどん刷って紙幣価値を下げ、
後者の道をまっしぐら。
「税金で借金を返せる難易度ランキング」
１位 日本 ２５７％
２位 スーダン ２１０％
３位 ギリシャ ２０７％
４位 エリトリア １７５％
５位 カーボベルデ １６１％
６位 イタリア １５５％
７位 スリナム １４１％
://twitter.com/fujimaki_takesi/status/1505469787174227973?s=20&t=3LeCpoO4fcZVopj2x8b9Dw
日本から始まる世界的株式市場の大暴落
それが最終的な暴落であることがはっきりするや否や、
マ仆レーヤは出現するでしょう。
マ仆レーヤが公に世界に現れるにつれて、
ＵＦＯがとてつもない数で姿を表すでしょう
(deleted an unsolicited ad)
989:１３２人目の素数さん
22/03/27 21:47:05.79 nv2FRWMk.net
質問いたします。よろしくお願いいたします。
ベクトル
a <- c(1, 2, 3, 3, 5, 5, 7) #1:7と連番にしたいが、一部値が重複している想定です
に対して、重複を修復するために
df <- data.frame(b=a, c=lag(a, default=a[1]-1))
apply(df, 1, function(x){ifelse(x[1]==x[2], x[1]+1, x[1])})
としましたが、applyを使っているのでかなり遅いです。
質問１：applyを使わない高速化できる手段はないでしょうか。
質問２：ifelseがベクトルに作用すると聞いたので（お恥ずかしいw）
ifelse(df$b == df$a,
としようとしたのですがエラーが出ました。
ifelseは異なる２つのベクトル（行数同じとしても）には作用できないのですね？
990:１３２人目の素数さん
22/03/27 21:55:41.37 C7ecChkw.net
>>941
何をどうしたいのか、具体的に書いてくれ。
何がしたいのか、さっぱりわからん。
c(1, 2, 3, 3, 5, 5, 7) を 1:7 にしたいだけなら、seq(along=a) でよいが、前と同じ値の場合は、前の値に1を足すというのであれば、rle を使えばなんとかなりそう。
991:１３２人目の素数さん
22/03/27 22:01:16.21 BmY4vlnd.net
質問が下手だと答えてもらえないぞ。
992:１３２人目の素数さん
22/03/27 22:03:58.44 nv2FRWMk.net
>>942
すいませんでした。詳しく書くと、
cベクトルにはcsvから読み込んだ、時間データの秒数が整数値で入っています。
(・csvデータには他に、時間データの年、月、日、時、分、秒、と
全て整数値がそれぞれ別の列で入っています。
・このcsvデータは１秒サンプリングです。
ただそのデータは秒がたまに重複しているのです。
重複している時のルールとして、直前の、１秒前の時間の秒が秒列に入っています。
例で挙げたように
1, 2, 3, 3, 5, 5, 7, 8, 9, 9
この重複を訂正したいのです。
そこで考えたのが 939で挙げたスクリプトだったのですが、
あまりにも遅く(1日分のデータ24*60*60行を訂正するのに）
質問いたしました。
993:１３２人目の素数さん
22/03/27 22:13:43.93 UiNoGcwm.net
おれもさっぱり質問の意図が理解できないが、重複除去なら集合演算の
b <- unique(a)
でbは
> b
[1] 1 2 3 5 7
になるよ。
994:939
22/03/27 22:17:25.42 nv2FRWMk.net
>>945
回答ありがとうございます。わかりづらくてすいません。
重複除去ではなく、重複のです。
942に書いたように、
秒データが
1,2,3,4,5,6,7,8
とあるべきところが
1,2,3,3,5,5,7,8
のように4秒のところが直前の３秒に、６秒のところが直前の５秒になっています。
これを修正したいのです。
995:１３２人目の素数さん
22/03/27 22:24:40.36 BmY4vlnd.net
そういうことか。
なら
a <- 1:8
で置換しちゃえばいいんじやないの？
てかまず国語を勉強しよう。
996:939
22/03/27 22:30:24.14 nv2FRWMk.net
>>947
ありがとうございます。
確かに国語力ないですね、今から考えると、、、
必要な前提を書ききれていません。
さらに必要な前提は、データが24*60*60秒分全部なく、１、２秒分抜けているのです。
997:１３２人目の素数さん
22/03/27 22:34:34.62 C7ecChkw.net
例が悪いね。
本当に1秒ごとに確実にデータが来ているなら、seq でも使って新たに値を作ればいいのだが、
きっと、途中でデータが来ない場合もあるのだろう。
例としては、c(1, 2, 2, 2, 9, 9, 11) のようなほうがよいのでは?
この場合、c(1, 2, 3, 4, 9, 10, 11) という答えになればよいのだろう。
r <- rle(a)
res <- lapply(seq(along=r$length), function(i) seq(r$values[i], len=r$length[i]))
do.call(c, res)
これでできると思うが、これで遅いなら、重複のところだけを修正するようにすればよい。
998:１３２人目の素数さん
22/03/27 22:35:09.61 UiNoGcwm.net
反省ゼロだなｗ
999:１３２人目の素数さん
22/03/27 22:37:05.89 BmY4vlnd.net
確認だけど、重複してる場合は必ず2つで、3つ以上はないの？
1000:939
22/03/27 22:38:26.76 nv2FRWMk.net
>>949
ありがとうございます。
参考に考えてみます。
>>950
返す返す申し訳ないです、、
1001:939
22/03/27 22:41:51.48 nv2FRWMk.net
>>951 データの特徴は、 １。重複は２つ。 ２。データの抜けは想定できない。 1日に５分ほどのこともあれば(付け加えて言うと、１日のデータをファイルに保存している) １日に1から３秒ほどの時もあります。(この時は、１日のデータをファイルに保存していない時)
1003:１３２人目の素数さん
22/03/27 22:42:10.78 UiNoGcwm.net
ヘタに抽象化せずに
隠さずに具体例を挙げて
FromとToを書けば無駄なやり取りをせずにすむ。
1004:１３２人目の素数さん
22/03/27 22:43:35.30 UiNoGcwm.net
>>953
そういうのが「ヘタな抽象化」です。
1005:１３２人目の素数さん
22/03/27 22:47:45.39 colv4Rae.net
まあ、落ち着け。時間データをクレンジングしたいんだとは思うけどもう少し元データについて整理しなされ。
1006:939
22/03/27 22:57:08.96 nv2FRWMk.net
>>955
>>956
ありがとうございます。
データが多いので、下手に抽象化してしまいましたが、
今手元のデータ例では、19:15:26のデータから秒列だけ上げると
下記のように19:15:36から秒列が重複することを始めます。
これが終了するのはデータによりけりで不明です。
26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 35, 37, 37, 39, 39, 41, 41, 43, 43, 45
1007:１３２人目の素数さん
22/03/27 22:58:55.54 UiNoGcwm.net
おれはもう脱落する。これは無理だ。
1008:１３２人目の素数さん
22/03/28 00:36:05.49 wfEbVN/P.net
これなら、修正箇所が少なくても多くてもパフォーマンスは変わらないはず。
comp.before <- function(x) {
# x の NA を直前の値で補完する.
valid <- !is.na(x) # NA でない(有効な)場所
v <- x[valid] # 有効な値
n <- diff(c(which(valid), length(x)+1)) # それを繰り返す回数
rep(v, n)
}
a <- c(1, 2, 2, 2, 9, 9, 11)
s <- seq(along=a) # 1から始まる等差数列
d <- ifelse(c(FALSE, diff(a) == 0), NA, s - a) # s から a の正しい部分を復元するための差
res <- s - comp.before(d) # 重複部分を修正するには、正しい部分の d を引き継げばよい
1009:１３２人目の素数さん
22/03/28 01:12:01.08 srY6abow.net
>>957
寝る前だったので、秒の部分だけでサンプルデータコードにしているので、そこは、
適宜、実際のデータに合わせて修正して。
library(tidyverse)
data.frame(
time = c(26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 35, 37, 37, 39, 39, 41, 41, 43, 43, 45)
) %>%
dplyr::mutate(diff = time - dplyr::lag(time, default = 0)) %>%
dplyr::mutate(new = dplyr::if_else(diff == 0, time + 1, time))
データ重複が2連続だけというのが前提です。3連続以上の場合は、最後の行の
処理を参考に追加してみてください。
1010:939
22/03/28 11:11:48.47 H+9Gp8uW.net
>>960
ありがとうございます。
なるほど、if_elseをそのように使うのですね。
条件の中で二つのベクトルを使ってダメだったので、諦めてました。
勉強になりました！！
1011:１３２人目の素数さん
22/03/28 11:56:48.44 srY6abow.net
>>961
lag()関数のdefault引数の指定があまりよくないので、
こちらの方が良いでしょう。処理手順を分かりやすく
するために処理を分割しています。
ただし、元データに欠損値（NA）がある場合は、意図した
結果を得られない場合がある点には注意してください。
data.frame(
time = c(26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 35, 37, 37, 39, 39, 41, 41, 43, 43, 45)
) %>%
dplyr::mutate(lag = dplyr::lag(time)) %>%
dplyr::mutate(diff = time - lag) %>%
dplyr::mutate(new = dplyr::case_when(diff == 0 ~ time + 1,
diff != 0 ~ time,
is.na(diff) ~ time))