18/04/26 18:58:06.22 UTg/BRA4.net
# URLﾘﾝｸ(www.igaku-shoin.co.jp)
# 相対リスク＝1.82を論文記載の対照群のリスク（36人/10万人）
# 1.82（95％信頼区間0.79-4.20）
rr=1.82
r0=36/10^5
r1=rr*r0
nn=seq(10^4,10^5,by=100)
rr2p.e <- function(n){
Epi::twoby2(matrix(c(n*r1,n*r0,n-n*r1,n-n*r0),ncol=2),print=FALSE)$p.value
}
plot(nn,sapply(nn,rr2p.e),type='l') ; abline(h=0.05, lty=3)
n=uniroot(function(n) rr2p.e(n)-0.05, c(10^4,10^5))$root
mat=matrix(c(n*r1,n*r0,n-n*r1,n-n*r0),ncol=2)
Epi::twoby2(mat)
491:卵の名無しさん
18/04/28 09:44:39.14 w/IREkFx.net
# 2 6 12 54 56+ 68 89 96 96 125+ 128+ 131+ 140+ 141+ 143 145+ 146 148+ 162+ 168 173+ 181+
time=c(2,6,12,54,56,68,89,96,96,125,128,131,140,141,143,145,146,148,162,168,173,181)
y=c(21/22,20/21,19/20,18/19,1,16/17,15/16,14/15,13/14,1,1,1,1,1,7/8,1,5/6,1,1,2/3,1,1)
survival=cumprod(y)
plot(time,survival,type='S',ylim=c(0,1))
censored=c(0,0,0,0,1,0,0,0,0,1,1,1,1,1,0,1,0,1,1,0,1,1)
event=1-censored
library(survival)
Surv(time,event)
plot(survfit(Surv(time,event)~1))
492:卵の名無しさん
18/04/28 11:58:57.59 w/IREkFx.net
カプランマイヤーのproduct limit法を理解していたらパッケージに頼らずにグラフも書けるな。
打ち切りポイントも描画できるように変更。
time=c(2,6,12,54,56,68,89,96,96,125,128,131,140,141,143,145,146,148,162,168,173,181)
y=c(21/22,20/21,19/20,18/19,1,16/17,15/16,14/15,13/14,1,1,1,1,1,7/8,1,5/6,1,1,2/3,1,1)
survival=cumprod(y)
censored=c(0,0,0,0,1,0,0,0,0,1,1,1,1,1,0,1,0,1,1,0,1,1)
plot(time,survival,type='s',ylim=c(0,1))
for(i in 1:22) if(censored[i]) points(time[i],survival [i],pch='+')
493:卵の名無しさん
18/04/30 11:29:32.46 T13xGV6f.net
# rule of thumb
# サンプルサイズと分散が同等な正規分布する集団からの無作為抽出で標本平均±標準誤差区間が重ならないときは母平均の有意差有無の判断はできない。
N=1000
f <- function(mx){
set.seed(123)
x1=rnorm(N,0 ,1)
x2=rnorm(N,mx,1)
t.test(x1,x2,var.equal = TRUE)$p.value
}
xx=seq(0,0.5,length=1000)
plot(xx,sapply(xx,f))
abline(h=0.05,lty=3)
uniroot(function(x) f(x)-0.05,c(0,1))$root
mx=0.0615 # 0.615以上で有意
g <- function(mx){
set.seed(123)
x1=rnorm(N,0 ,1)
x2=rnorm(N,mx,1)
n1=length(x1) ; n2=length(x2)
m1=mean(x1) ; m2=mean(x2)
s1=sd(x1) ; s2=sd(x2)
se1=s1/sqrt(n1) ; se2=s2/sqrt(n1)
(m1+se1) > (m2-se2) # overlap?
}
xx=seq(0,0.5,length=1000)
plot(xx,sapply(xx,g))
abline(h=0,lty=3)
uniroot(g,c(0,1))$root # non-overlap には0.0369以上
f(0.07) # 有意差あり
g(0.07) # overlapなし
f(0.05) # 有意差なし
g(0.05) # overlapなし
494:卵の名無しさん
18/04/30 12:22:04.07 T13xGV6f.net
# rule of thumb
# サンプルサイズと分散が同等な正規分布する集団からの無作為抽出で
# 標本平均±標準誤差区間が重なるときは母平均に有意差はない。
##
495:卵の名無しさん
18/04/30 14:24:19.80 T13xGV6f.net
# rule of thumb
# サンプルサイズと分散が同等な正規分布する集団からの無作為抽出で
# 標本平均±標準誤差区間が重なるときは母平均に有意差はない(p>0.05)。
#
x1=rnorm(n,m1,s)
x2=rnorm(n,m2,s) # m1 < m2
t.stat=(m1-m2)/(sqrt(2/n)*s) # < 0
=(m1-m2)*sqrt(n)/(sqrt(2)*s)
m1 + se > m2 -se # overlap
m1 - m2 > -2*se = -2*s/sqrt(n)
m1 - m2 = t.stat*(sqrt(2)*s)/sqrt(n) > -2*s/sqrt(n)
t.stat > -2/sqrt(2)=-sqrt(2)
df=198 # 例
curve(pt(x,df),-5,0) # x < 0で増加関数
pt(t.stat,df) > pt(-sqrt(2),df)
pt(-sqrt(2),df)
# dfを変化させる
curve(pt(-sqrt(2),x),200,ylim=c(0,0.5),xlab='df',ylab='p.value',lwd=2)
abline(h=0.05,lty=3)
pt(-sqrt(2),Inf) ; pnorm(-sqrt(2))
text(0,0.05,'0.05')
496:卵の名無しさん
18/04/30 14:24:33.69 T13xGV6f.net
# 平均値の95%信頼区間が重なった場合にはそれで有意差判定はできない。
set.seed(123)
x1=rnorm(100,0 ,1)
x2=rnorm(100,0.3 , 1)
t.test(x1,x2,var.equal = TRUE)$p.value # > 0.05
t.test(x1)$conf[2] > t.test(x2)$conf[1] # overlap TRUE
plot(NULL,NULL,ylim=c(0,0.05),xlim=c(-0.75,0.75),yaxt='n',ann=FALSE,bty='l')
segments(t.test(x1)$conf[1],0,t.test(x1)$conf[2],0,lwd=5,col=1)
segments(t.test(x2)$conf[1],0.01,t.test(x2)$conf[2],0.01,lwd=5,col=2)
set.seed(123)
x1=rnorm(100,0 ,1)
x2=rnorm(100,0.5 ,1)
t.test(x1,x2,var.equal = TRUE)$p.value # < 0.05
t.test(x1)$conf[2] > t.test(x2)$conf[1] # overlap TRUE
plot(NULL,NULL,ylim=c(0,0.05),xlim=c(-0.75,0.75),yaxt='n',ann=FALSE,bty='l')
segments(t.test(x1)$conf[1],0,t.test(x1)$conf[2],0,lwd=5,col=1)
segments(t.test(x2)$conf[1],0.01,t.test(x2)$conf[2],0.01,lwd=5,col=2)
497:卵の名無しさん
18/05/01 08:02:38.81 v/zhryAm.net
T.test2=function(n,dm,var1,var2){
SE12=sqrt((1/n+1/n)*((n-1)*var1+(n-1)*var2)/((n-1)+(n-1)))
T=dm/SE12
2*pt(abs(T),n-1+n-1,lower.tail = FALSE)
}
m1-m2=dm
n1=n2=n
T.test3=function(n,dm,var1,var2){
SE12=sqrt((2/n)*(var1+var2)/2))
T=dm/SE12
2*pt(abs(T),n-1+n-1,lower.tail = FALSE)
}
T.test4=function(x,n=1000,dm=1,var1=1){
SE=sqrt((2/n)*(var1+x)/2)
T=dm/SE
2*pt(abs(T),n-1+n-1,lower.tail = FALSE)
}
curve(T.test4(x),0,100,xlab='variance=x against variance=1',ylab='p.value')
498:卵の名無しさん
18/05/01 08:23:44.83 v/zhryAm.net
# t検定(生データなし,等分散不問)
Welch.test=function(n1,n2,m1,m2,var1,var2){
T=(m1-m2)/sqrt(var1/n1+var2/n2)
df=(var1/n1+var2/n2)^2 / (var1^2/n1^2/(n1-1)+var2^2/n2^2/(n2-1))
p.value=2*pt(abs(T),df,lower.tail = FALSE)
return(p.value)
}
n1=n2=n
m1-m2=dm
var2=x
Welch.test2=function(x,n=1000,dm=1,var1=1){
T=dm/sqrt(var1/n+x/n)
df=(var1/n+x/n)^2 / (var1^2/n^2/(n-1)+x^2/n^2/(n-1))
p.value=2*pt(abs(T),df,lower.tail = FALSE)
return(p.value)
}
curve(Welch.test2(x),0,100,xlab='variance=x against variance=1',ylab='p.value')
499:卵の名無しさん
18/05/01 09:37:21.50 UP3hBRO4.net
N=1000
n1=40
n2=60
SDR=10
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (10^4,f2.b())
hist (p.t)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (10^4,f2.bW())
hist (p.W)
500:卵の名無しさん
18/05/01 10:15:19.84 UP3hBRO4.net
N=1000
n1=40
n2=60
SDR=1
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (10^4,f2.b())
hist (p.t)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (10^4,f2.bW())
hist (p.W)
501:卵の名無しさん
18/05/01 10:21:41.36 UP3hBRO4.net
N=1000
n1=40
n2=60
SDR=1
dm=0.1
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR+dm
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (10^4,f2.b())
hist (p.t)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (10^4,f2.bW())
hist (p.W)
502:卵の名無しさん
18/05/01 12:51:04.77 UP3hBRO4.net
T.test4=function(var2,n1,n2,dm=0.5,var1=1){
SE12=sqrt((1/n1+1/n2)*((n1-1)*var1+(n2-1)*var2)/((n1-1)+(n2-1)))
T=dm/SE12
2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
curve(T.test4(x,40,160,0.1),0,100,xlab='variance=x against variance=1',ylab='p.value')
503:卵の名無しさん
18/05/01 14:51:52.52 W1Ln0POm.net
# URLﾘﾝｸ(www.rips-irsp.com)
# Why Psychologists Should by Default Use Welch’s t-test
# Instead of Student’s t-test
par(mfrow=c(2,1))
N=1000
n1=10
n2=90
SDR=10
dm=1
k=10^4
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR+dm
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (k,f2.b())
hist (p.t,main='Student\'s t.test',col=sample(colours(),1))
mean(p.t < 0.05) # power(when dm!=0) or Type I error(when dm=0)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (k,f2.bW())
hist (p.W,main='Welch\'s t.test',col=sample(colours(),1))
mean(p.W < 0.05) # power(when dm!=0) or Type I error(when dm=0)
504:卵の名無しさん
18/05/01 14:57:05.14 W1Ln0POm.net
Why Psychologists Should by Default Use Welch’s t-test Instead of Student’s t-test
URLﾘﾝｸ(www.rips-irsp.com)
に触発されてシミュレーションの条件（サンプルサイズや分散比）を変えて追試してみた。
結果にびっくり
URLﾘﾝｸ(i.imgur.com)
あたかもド底辺シリツ医大卒を最高学府を履修したと呼ぶほどの差異に等しい。
Welchのrobustnessに再度感銘した。
505:卵の名無しさん
18/05/01 19:34:12.07 W1Ln0POm.net
par(mfrow=c(2,1))
N=1000
n1=50
n2=25
SDR=5
dm=0 # null hypothesis is true when dm==0
k=10^4
set.seed(123)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR+dm
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (k,f2.b())
hist (p.t,main='Student\'s t.test',col=sample(colours(),1),
xlab=paste('n1 = ',n1,', sd1 = 1 ; n2 = ',n2,' , sd2 = ',SDR))
mean(p.t < 0.05) # power(when dm!=0) or Type I error(when dm=0)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (k,f2.bW())
hist (p.W,main='Welch\'s t.test',col=sample(colours(),1),
xlab=paste('n1 = ',n1,', sd1 = 1 ; n2 = ',n2,' , sd2 = ',SDR))
mean(p.W < 0.05) # power(when dm!=0) or Type I error(when dm=0)
506:卵の名無しさん
18/05/04 16:46:09.21 6BNcb2/S.net
x=log.nor=seq(-8.0,-3.0,by=0.5) ; n=length(x)
y=relax=c(2.6,10.5,15.8,21.1,36.8,57.9,73.7,89.5,94.7,100.0,100.0)
(re.nls=nls(y ~ Top/(1+10^((LogEC50-x)*HillSlope)),
start=c(Top=100,LogEC50=-5,HillSlope=1),algo='port'))
nls2R2(re.nls,y)
plot(x,y,bty='l',pch=19)
lines(x,predict(re.nls))
#
R2 <- function(dat){
n = nrow(dat)
sm = summary(
nls(dat$y ~ Top/(1+10^((LogEC50-dat$x)*HillSlope)),
start=c(Top=100,LogEC50=-5,HillSlope=0.5),algorithm = 'port')
)
SSalt = sum(sm$residuals^2)
SSnull = var(dat$y)*(n-1)
Rsq=1-SSalt/SSnull
Rsq
# k=sm$df[1]
# AdjRsq=1-(1-Rsq)*(n-1)/(n-k-1)
# data.frame(Rsq,AdjRsq)
}
xy
507:=data.frame(x,y) R2(xy) library(boot) re=boot(xy,function(data,indices) R2(data[indices,]),R=1000) re$t0 BEST::plotPost(re$t) ; HDInterval::hdi(re$t) quantile(re$t,c(.025,.975)) nls2R2(re.nls,y)
508:卵の名無しさん
18/05/04 17:55:57.07 6BNcb2/S.net
anova(lm)$Pr[1])
min( pf(summary(lm)$fstatistic[1],
summary(lm)$fstatistic[2],summary(lm)$fstatistic[2],lower=FALSE),
pf(summary(lm)$fstatistic[1],
summary(lm)$fstatistic[2],summary(lm)$fstatistic[2] )
509:卵の名無しさん
18/05/04 18:01:55.70 6BNcb2/S.net
lm = lm(formula,data)
lm2p <- function(lm){ # p_value=anova(lm)$Pr[1]
f=summary(lm)$fstatistic
min(
pf(f[1],f[2l,f[3],lower=FALSE),
pf(f[1],f[2l,f[3])
)
}
510:卵の名無しさん
18/05/04 20:28:04.76 6BNcb2/S.net
＃＃ ｈｔｔｐ：／／ｗｗｗ．ｐｕｂｌｉｃ．ａｓｕ．ｅｄｕ／～ｇａｓｗｅｅｔｅ／ｃｒｊ６０４／ｒｅａｄｉｎｇｓ／１９８３－Ｆｒｅｅｄｍａｎ％２０（Ｓｃｒｅｅｎｉｎｇ％２０Ｒｅｇｒｅｓｓｉｏｎ％２０Ｅｑｕａｔｉｏｎｓ）．ｐｄｆ
ｎｏｉｓｅ２ｓｉｇ ＜－ ｆｕｎｃｔｉｏｎ（ｓｅｅｄ，ｓｅｑｕｅｎｔｉａｌ＝ＦＡＬＳＥ）｛
ｓｅｔ．ｓｅｅｄ（ｓｅｅｄ）
Ｎｏｉｓｅ＝ｍａｔｒｉｘ（ｒｎｏｒｍ（１００＊５１），ｎｃｏｌ＝５１，ｎｒｏｗ＝１００）
＃ ｆｉｒｓｔ ｐａｓｓ
ｌｍ１＝ｌｍ（Ｎｏｉｓｅ［，５１］ ～ Ｎｏｉｓｅ［，－５１］） ＃ Ｎｏｉｓｅ［，５１］ ： 目的変数Ｙ
ｃａｔ（＇\ｎＲ２ （１ｓｔ） ＝ ＇，ｓｕｍｍａｒｙ（ｌｍ１）＄ｒ．ｓｑｕａｒｅｄ）
ｃａｔ（＇\ｎｐ．ｖａｌｕｅ （１ｓｔ） ＝ ＇，ａｎｏｖａ（ｌｍ１）＄Ｐｒ［１］）
ｃｏｅｆｓ１＝ｓｕｍｍａｒｙ（ｌｍ１）＄ｃｏｅｆ［，＇Ｐｒ（＞｜ｔ｜）＇］
ｌｅｎｇｔｈ（ｃｏｅｆｓ１）
ｃｏｅｆｓ１［１］
ｃｏｅｆｆｓ１＝ｃｏｅｆｓ１［－１］
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．２５ ＝＇，ｓｕｍ（ｃｏｅｆｆｓ１＜０．２５））
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．０５ ＝＇，ｓｕｍ（ｃｏｅｆｆｓ１＜０．０５））
ｉｎｄｘ２５＝ｗｈｉｃｈ（ｃｏｅｆｆｓ１＜０．２５）
511:卵の名無しさん
18/05/04 20:28:19.24 6BNcb2/S.net
＃ ｓｅｃｏｎｄ ｐａｓｓ
ｌｍ２＝ｌｍ（Ｎｏｉｓｅ［，５１］ ～ Ｎｏｉｓｅ［，ｉｎｄｘ２５］）
ｃａｔ（＇\ｎ\ｎＲ２ （２ｎｄ） ＝ ＇，ｓｕｍｍａｒｙ（ｌｍ２）＄ｒ．ｓｑｕａｒｅｄ）
ｃａｔ（＇\ｎｐ．ｖａｌｕｅ（２ｎｄ） ＝ ＇，ａｎｏｖａ（ｌｍ２）＄Ｐｒ［１］）
ｃｏｅｆｓ２＝ｓｕｍｍａｒｙ（ｌｍ２）＄ｃｏｅｆ［，＇Ｐｒ（＞｜ｔ｜）＇］
ｌｅｎｇｔｈ（ｃｏｅｆｓ２）
ｃｏｅｆｓ２［１］
ｃｏｅｆｆｓ２＝ｃｏｅｆｓ２［－１］
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．２５ （２ｎｄ） ＝＇，ｓｕｍ（ｃｏｅｆｆｓ２＜０．２５））
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．０５ （２ｎｄ） ＝＇，ｓｕｍ（ｃｏｅｆｆｓ２＜０．０５））
ｃａｔ（＇\ｎ\ｎ＇）
ｉｆ（ｓｅｑｕｅｎｔｉａｌ）｛
ｃａｔ（＇Ｈｉｔ Ｒｅｔｕｒｎ Ｋｅｙ ｉｎ ｃｏｎｓｏｌｅ ｗｉｎｄｏｗ＇）
ｎｏ＿ｓａｖｅ ＜－ ｓｃａｎ（ｎ＝１， ｗｈａｔ＝ｃｈａｒａｃｔｅｒ（）， ｑｕｉｅｔ＝ＴＲＵＥ）
｝
｝
ｎｏｉｓｅ２ｓｉｇ（１）
ｆｏｒ（ｉ ｉｎ ２：５） ｎｏｉｓｅ２ｓｉｇ（ｉ，ｓｅｑ＝ＴＲＵＥ）
512:卵の名無しさん
18/05/07 08:41:59.73 /y3oEDZQ.net
UG.Purge <- function(alt='two.sided'){
pv=numeric(10)
for(j in 1:10){
L=list()
tmp=UG[-j,]
for(i in 1:6) L=append(L,list(tmp[,i]))
pv[j]=jonckheere(L,cat=FALSE,alternative = alt)
}
return(pv)
}
which(UG.Purge('two') < 0.05)
which(UG.Purge('increasing') < 0.05)
513:卵の名無しさん
18/05/08 17:13:30.89 4Ru4wKk7.net
＃ ｈｔｔｐ：／／ｆｉｌｅ．ｓｃｒａｔｃｈｈｉｔ．ｐａｚｒｕ．ｃｏｍ／ＴｋｙＳｃｈＢｎｆ．ｔｘｔ
＃＃ ｆｉｔｔｅｄ ｉｓ ａ ｓｅｑｕｅｎｃｅ ｏｆ ｍｅａｎｓ
＃＃ ｎｉｓ ｉｓ ａ ｃｏｒｒｅｓｐｏｎｄｉｎｇ ｓｅｑｕｅｎｃｅ ｏｆ ｓａｍｐｌｅ ｓｉｚｅｓ ｆｏｒ ｅａｃｈ ｍｅａｎ
＃＃ ｄｆ ｉｓ ｔｈｅ ｒｅｓｉｄｕａｌ ｄｆ ｆｒｏｍ ｔｈｅ ＡＮＯＶＡ ｔａｂｌｅ
＃＃ ＭＳＥ ＝ ｍｅａｎ ｓｑｕａｒｅｄ ｅｒｒｏｒ ｆｒｏｍ ｔｈｅ ＡＮＯＶＡ ｔａｂｌｅ
＃＃ ｃｏｎｆ．ｌｅｖｅｌ ｉｓ ｔｈｅ ｆａｍｉｌｙ－ｗｉｓｅ ｃｏｎｆｉｄｅｎｃｅ ｌｅｖｅｌ， ｄｅｆａｕｌｔｓ ｔｏ ．９５
ｐａｉｒｗｉｓｅ．ｓｃｈｅｆｆｅ．ｔ�
514:�ｓｔ ＜－ ｆｕｎｃｔｉｏｎ（ｒｅ．ａｏｖ， ｃｏｎｆ．ｌｅｖｅｌ ＝ ０．９５）｛ ｍｏｄｅｌ ＝ ｒｅ．ａｏｖ＄ｍｏｄｅｌ ｃｏｌｎａｍｅｓ（ｍｏｄｅｌ） ＝ ｃ（＇ｓｃｏｒｅ＇，＇ｇｒｏｕｐ＇） ｆｉｔｔｅｄ ＝ ｗｉｔｈ（ｍｏｄｅｌ，ｔａｐｐｌｙ（ｓｃｏｒｅ，ｇｒｏｕｐ，ｍｅａｎ）） ｎｉｓ ＝ ｗｉｔｈ（ｍｏｄｅｌ，ｔａｐｐｌｙ（ｓｃｏｒｅ，ｇｒｏｕｐ，ｌｅｎｇｔｈ）） ｄｆ ＝ ｒｅ．ａｏｖ＄ｄｆ．ｒｅｓｉｄｕａｌ ＭＳＥ ＝ ｓｕｍｍａｒｙ（ｒｅ．ａｏｖ）［［１］］［＇Ｒｅｓｉｄｕａｌｓ＇，＇Ｍｅａｎ Ｓｑ＇］ ｒ ＝ ｌｅｎｇｔｈ（ｆｉｔｔｅｄ） ｐａｉｒｓ ＝ ｃｏｍｂｎ（１：ｒ，２） ｄｉｆｆｓ ＝ ｆｉｔｔｅｄ［ｐａｉｒｓ［１，］］ － ｆｉｔｔｅｄ［ｐａｉｒｓ［２，］］
515:卵の名無しさん
18/05/08 17:13:57.38 4Ru4wKk7.net
Ｔ ＝ ｓｑｒｔ（（ｒ－１）＊ｑｆ（ｃｏｎｆ．ｌｅｖｅｌ，ｒ－１，ｄｆ））
ｈｗｉｄｔｈｓ ＝ Ｔ＊ｓｑｒｔ（ＭＳＥ＊（１／ｎｉｓ［ｐａｉｒｓ［１，］］ ＋ １／ｎｉｓ［ｐａｉｒｓ［２，］］））
ｆｖｓ ＝ （ｄｉｆｆｓ＾２）／（ＭＳＥ＊（１／ｎｉｓ［ｐａｉｒｓ［１，］］ ＋ １／ｎｉｓ［ｐａｉｒｓ［２，］］））／（ｒ－１）
ｐｖｓ ＝ １－ｐｆ（ｆｖｓ， ｒ－１， ｄｆ）
ｖａｌ ＝ ｃｂｉｎｄ（ｄｉｆｆｓ － ｈｗｉｄｔｈｓ， ｄｉｆｆｓ， ｄｉｆｆｓ ＋ ｈｗｉｄｔｈｓ， ｆｖｓ， ｒ－１， ｄｆ， ｐｖｓ）
ｄｉｍｎａｍｅｓ（ｖａｌ） ＝ ｌｉｓｔ（ｐａｓｔｅ（＂ｓｕｂｊｅｃｔ＂，ｐａｉｒｓ［１，］，＂ － ｓｕｂｊｅｃｔ＂， ｐａｉｒｓ［２，］，
ｓｅｐ ＝ ＂＂）， ｃ（＂Ｌｏｗｅｒ＂， ＂Ｄｉｆｆ＂，＂Ｕｐｐｅｒ＂， ＂Ｆ Ｖａｌｕｅ＂， ＂ｎｍｄｆ＂， ＂ｄｎｄｆ＂， ＂Ｐｒ（＞Ｆ）＂））
Ｓｃ ＝ ａｓ．ｍａｔｒｉｘ（ｖａｌ）
ｐ．Ｓｃ ＝ Ｓｃ［，＇Ｐｒ（＞Ｆ）＇］
ｎ ＝ ｒｅ．ａｏｖ＄ｒａｎｋ － １
ｍａｔ．Ｓｃ ＝ ｍａｔｒｉｘ（ｒｅｐ（ＮＡ，ｎ＊ｎ），ｎ）
ｆｏｒ（ｋ ｉｎ １：ｎ） ｍａｔ．Ｓｃ［，ｋ］ ＝ ｃ（ｒｅｐ（ＮＡ，ｋ－１）， ｐ．Ｓｃ［（（ｋ－１）＊（ｎ＋１－ｋ／２）＋１）：（（ｋ－１）＊（ｎ＋１－ｋ／２）＋１＋ｎ－ｋ）］）
ｃｏｌｎａｍｅｓ（ｍａｔ．Ｓｃ） ＝ ｍｏｄｅｌ［［２］］［１：ｎ］
ｒｏｗｎａｍｅｓ（ｍａｔ．Ｓｃ） ＝ ｍｏｄｅｌ［［２］］［２：（ｎ＋１）］
ｄ ＝ ｒｏｕｎｄ（ｍａｔ．Ｓｃ，５）
ｄ［ｉｓ．ｎａ（ｄ）］ ＝ ＇－＇
ｐｒｉｎｔ（ｄ，ｑｕｏｔｅ ＝ ＦＡＬＳＥ）
ｉｎｖｉｓｉｂｌｅ（ｖａｌ）
｝
516:卵の名無しさん
18/05/09 14:26:52.13 L5gMx8M1.net
# disese, No disease
# positive test TP(A) FP(B)
# negative test FN(C) TN(D)
PV <- function(prevalence=c(10^-4,1),sn=0.80,sp=0.95,...){
plot(NULL,NULL,xlim=prevalence,ylim=c(0,1),xlab='Prevalence',
ylab='Predicative Value(Test Credibility)',type='n',bty='l',
main=paste0('感度 = ',sn, ' 特異度 = ',sp), ...)
legend('center',bty='n',lwd=2,lty=1:3,legend=c('陽性適中率','陰性適中率','偽陽性率'),
col=c('red','darkgreen','brown'))
ppv<-function(prevalence,sensitivity=sn,specificity=sp){
PPV=prevalence*sensitivity/
(prevalence*sensitivity+(1-prevalence)*(1-specificity))
return(PPV)
}
curve(ppv(x),lty=1,lwd=2,col='red', add=TRUE)
npv<-function(prevalence,sensitivity=sn,specificity=sp){
NPV=(1-prevalence)*specificity/
( (1-prevalence)*specificity + prevalence*(1-sensitivity) )
return(NPV)
}
curve(npv(x),lty=2,lwd=2,col='darkgreen',add=TRUE)
false_negative_rate <- function(prevalence,sensitivity=sn,specificity=sp){
FNR <- prevalence*(1-sensitivity)/
( (1-prevalence)*specificity + prevalence*(1-sensitivity) )
return(FNR)
}
curve(false_negative_rate(x), lty=3,lwd=2,col='brown',add=TRUE)
}
517:卵の名無しさん
18/05/11 20:34:40.05 MdgPLO2C.net
zodiac=c('Aquarius','Aries','Cancer','Capricorn','Gemini','Leo','Libra','Pisces','Sagittarius','Scorpio','Taurus','Virgo')
N=c(856301,888348,917553,844635,937615,903009,897503,893332,846813,850128,918512,921196)
X=c(1433,1476,1496,1343,1553,1497,1350,1522,1277,1297,1534,1445)
iP=which(zodiac=='Pisces')
p.Pisces=prop.test(c(X[iP],sum(X[-iP])),c(N[iP],sum(N[-iP])))$p.value
p0=sum(X)/sum(N)
se
518:t.seed(16062006) p.max=numeric() for(i in 1:10^4){ xi=rbinom(length(N),N,p0) I=which.max(xi/N) re=prop.test(c(xi[I],sum(xi[-I])),c(N[I],sum(N[-I]))) p.max[i]=re$p.value } mean(p.max <= p.Pisces) cat("P-value for Pisces and CHF: ", p.max, file="Pisces.txt",fill=TRUE, append=TRUE)
519:卵の名無しさん
18/05/13 14:56:25.30 lYhicRHD.net
TRAP9p <- function(seed=NULL){
set.seed(seed)
A0=rnorm(10,50,10)
A1=A0+rnorm(10,0,5)+5
B0=rnorm(10,45,10)
B1=B0+rnorm(10,0,5)+5
A0=round(A0)
B0=round(B0)
A1=round(A1)
B1=round(B1)
a=t.test(A0,A1,paired=TRUE)$p. # 外部有意差なし
b=t.test(B0,B1,paired=TRUE)$p. # 内部有意差あり
ab=t.test(A1-A0,B1-B0)$p. # 加点差有意差なし
a0b0=t.test(A0,B0)$p.# base差なし
a1b1=t.test(A1,B1)$p.# final差
list(p=c(外部=a,内部=b,差異=ab,base差=a0b0,final差=a1b1),A0=A0,B0=B0,A1=A1,B1=B1)
}
TRAP9p()$p
rep=replicate(1000,TRAP9p()$p)
mean(rep['外部',]>0.05 & rep['内部',]<0.05 & rep['差異',]>0.05
& re['base差',] > 0.05 & re['final差',]<0.05)
seeds=c(4799, 10638, 12173, 12908, 13671, 17145, 18955)
re=sapply(seeds,function(x) TRAP9p(x)$p) # wait several ten seconds
idx=which(re['外部',]>0.05 & re['内部',]<0.05 & re['差異',]>0.05
& re['base差',] > 0.05 & re['final差',]<0.05)
idx ; re[,idx]
Tx=as.data.frame(TRAP9p(4799)[2:5]) # 4799 10638 12173 12908 13671 17145 18955
U=with(Tx,data.frame(外部.無=A0,外部.有=A1,内部.無=B0,内部.有=B1))
520:卵の名無しさん
18/05/13 14:57:15.49 lYhicRHD.net
TRAP9p()$p
rep=replicate(1000,TRAP9p()$p)
mean(rep['外部',]>0.05 & rep['内部',]<0.05 & rep['差異',]>0.05
& re['base差',] > 0.05 & re['final差',]<0.05)
seeds=c(4799, 10638, 12173, 12908, 13671, 17145, 18955)
re=sapply(seeds,function(x) TRAP9p(x)$p) # wait several ten seconds
idx=which(re['外部',]>0.05 & re['内部',]<0.05 & re['差異',]>0.05
& re['base差',] > 0.05 & re['final差',]<0.05)
idx ; re[,idx]
Tx=as.data.frame(TRAP9p(4799)[2:5]) # 4799 10638 12173 12908 13671 17145 18955
U=with(Tx,data.frame(外部.無=A0,外部.有=A1,内部.無=B0,内部.有=B1))
U
.plot = function(){
plot (NULL,xlim=c(0,5),ylim=c(30,80),bty='l',type='n',
xaxt='n',ylab='score',xlab='',main='「任意」の寄付金効果')
axis(side=1, at=c(1.5,3.5), labels=c('外部入学','内部進学'))
points(rep(1,10),U[,1])
points(rep(2,10),U[,2])
segments(rep(1,10),U[,1],rep(2,10),U[,2])
points(rep(3,10),U[,3],pch=19)
points(rep(4,10),U[,4],pch=19)
segments(rep(3,10),U[,3],rep(4,10),U[,4])
} ; .plot()
with(U,t.test(外部.無,外部.有,paired=TRUE))$p.value # 外部での効果
with(U,t.test(内部.無,内部.有,paired=TRUE))$p.value # 内部での効果
with(U,t.test(外部.有-外部.無,内部.有-内部.無)$p.value) # 内外部効果比較
with(U,t.test(log(外部.有/外部.無),log(内部.有/内部.無))$p.value)
521:卵の名無しさん
18/05/14 08:43:54.12 R8CQudRE.net
alpha=0.05 # Pr(significantH0)
power=0.80 # Pr(significant|H1)
N=100
prior=0.10
disease=prior*N
TP=disease*power
FP=N*(1-prior)*alpha
(FPRP=FP/(TP+FP))
alpha.L=0.045
alpha.U=0.055
522:卵の名無しさん
18/05/14 14:11:13.07 elEQ53K0.net
sensitivity=0.80 # Pr(reject|H1)
specificity=1-0.05
alpha=1-specificity. # Pr(reject|H0)
n=1000
prior=0.1
prop.p<- function (){
a=rbinom(1,n, prior)
b=rbinom(1,n, prior)
prop.test(c(a,b),c(n,n))$p.value
}
k=1000
re=replicate (k,prop.p()1)
hist(re)
mean(re<0.05)
mean (0.045<re & re<0.050)
TP=sensitivfity*prior
do qFP=(1-prior)*(1-sensitivity)
FPRP=FP/(TP+FP)
523:卵の名無しさん
18/05/15 08:47:26.31 NGCOKTk3.net
knowing that the data are ‘rare’ when there is no true difference is of little use unless one determines whether
or not they are also ‘rare’ when there is a true difference.
524:卵の名無しさん
18/05/15 18:26:03.26 NGCOKTk3.net
BayesFactor=Pr(sig|H0)/Pr(sig|H1)=alpha/power=(1-specificity)/sensitivity=FP/TP=1/pLR
pLR=TP/FP=sensitivity/(1-specificity)=1/BayesFactor
nLR=FN/TN=(1-sensitivity)/specificity
525:卵の名無しさん
18/05/16 07:16:23.92 nUz5W2k6.net
calc.FPR.p <- function(r1,r2,n1,n2,alpha=0.05){ #
526:n1=n2 p.val=prop.test(c(r1,r2),c(n1,n2))$p.value k=2 df=k-1 Pi=c(r1/n1,r2/n2) n=mean(n1,n2) theta.i=asin(sqrt(Pi)) delta=4*var(theta.i)*(k-1) ncp=n*delta power=pchisq(qchisq(1-alpha,df),df,ncp,lower=FALSE) qcrit=qchisq(max(p.val,1-p.val),df,ncp=0) curve(dchisq(x,df),0,15,ann=FALSE,bty='n') # H0 curve(dchisq(x,df,ncp),add=TRUE, lty=2) # H1 abline(v=qcrit,lty=3) x0=qcrit y0=dt(x0,df,0) x1=x0 y1=dchisq(x1,df,ncp=ncp) FPR=y0/(y0+y1) print(c(p.value=p.val,FPR=FPR),digits=3) } calc.FPR.p(85,95,100,100)
527:卵の名無しさん
18/05/16 20:15:48.59 nUz5W2k6.net
calc.FPR.chisq <- function(r,alpha=0.05){ # r=c(37,21,21,21) vs n=c(25,25,25,25)
m=length(r)
n=rep(mean(r),m)
p.val=chisq.test(rbind(r,n))$p.value
df=m-1
P0=n/sum(n)
P1=(r/n)/sum(r/n)
N=sum(r) # == sum(n)
ncp1=0
for(i in 1:m) ncp1=ncp1+N*(P1[i]-P0[i])^2/P0[i]
qcrit=qchisq(1-alpha,df,ncp=0)
curve(dchisq(x,df),0,20,xlab=quote(chi),ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
abline(v=qcrit,col='gray') ; text(qcrit,0,round(qcrit,2))
legend('topright',bty='n',legend=c('H0:Central','H1:Noncentral'),col=1:2,lty=1:2)
power=pchisq(qcrit,df,ncp1,lower=FALSE)
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR=y0/(y0+y1)
print(c(power=power,p.value=p.val,FPR=FPR),digits=3)
}
calc.FPR.chisq(c(10,16,34,9,10,26))
calc.FPR.chisq(c(1,0,9,2,1,1))
528:卵の名無しさん
18/05/17 06:32:52.54 1zxWbYbO.net
サンプルサイズが異なっても算出できるように改造した。
パーッケージpwrのpwr.2p2n.testの.コードを参考にした。
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(r1/n1,r2/n2)
theta.i=asin(sqrt(Pi))
delta=4*var(theta.i)*(k-1)
N=2/(1/n1+1/n2)
ncp1=N*delta
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
curve(dchisq(x,df),0,20,xlab='ChiSquare',ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
abline(v=qcrit,col='gray')
legend('top',bty='n',legend=c('H0','H1'),col=1:2,lty=1:2)
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR=y0/(y0+y1)
VAL=c(power=power.chi,p.value=p.val,FPR=FPR)
print(VAL,digits=3)
invisible(VAL)
}
529:卵の名無しさん
18/05/17 07:09:30.35 1zxWbYbO.net
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(r1/n1,r2/n2)
theta.i=asin(sqrt(Pi))
delta=4*var(theta.i)*(k-1)
N=2/(1/n1+1/n2)
ncp1=N*delta
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
curve(dchisq(x,df),0,ncp1*3,xlab=quote(chi),ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
abline(v=qcrit,col='gray')
legend('top',bty='n',legend=c('H0','H1'),col=1:2,lty=1:2)
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR.equal=y0/(y0+y1)
FPR.less=p.val/(p.val+power.chi)
VAL=c(power=power.chi,p.value=p.val,FPR.equal=FPR.equal,FPR.less=FPR.less)
print(VAL,digits=3)
invisible(VAL)
}
calc.FPR.p2(95,85,100,110)
530:卵の名無しさん
18/05/17 10:16:58.64 Z8D8umCF.net
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(p1,p2)
theta.i=asin(sqrt(Pi)) # arcsine conversion
delta=4*var(theta.i)*(k-1) # sum of squares
N=2/(1/n1+1/n2) # harmonic mean, subcontrary mean
ncp1=N*delta # non-central parameter
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
# qchisq(1-0.05,1): 3.841
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
# qcrit = prop.test(c(r1,r2),c(n1,n2))$statistic
curve(dchisq(x,df),0,2*qcrit,xlab=quote(chi),ylab='Density',bty='n') # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2) # H1
# power.chi=AUC of right half of the H1 curve
abline(v=qcrit,col='gray')
legend('top',bty='n',legend=c('H0','H1'),col=1:2,lty=1:2)
text(qcrit,0,round(qcrit,2))
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR.equal=y0/(y0+y1) # length ratio
FPR.less=p.val/(p.val+power.chi) # area ratio
FPR.alpha=alpha/(alpha+power.chi) # FPR before analysis
VAL=c(power=power.chi,p.value=p.val,FPR.equal=FPR.equal,
FPR.less=FPR.less,FPR.alpha=FPR.alpha)
print(VAL,digits=3)
invisible(VAL)
}
calc.FPR.p2(85,95,100,100, alpha=0.05)
calc.FPR.p2(85,95,100,100, alpha=0.01)
531:卵の名無しさん
18/05/17 11:17:41.13 Z8D8umCF.net
calc.FPR.p2 <- function(r1,r2,n1,n2,alpha=0.05){
p.val=prop.test(c(r1,r2),c(n1,n2))$p.value
k=length(c(r1,r2))
df=k-1
p1=r1/n1 ; p2=r2/n2
Pi=c(p1,p2)
theta.i=asin(sqrt(Pi)) # arcsine conversion
delta=4*var(theta.i)*(k-1) # sum of squares
N=2/(1/n1+1/n2) # harmonic mean, subcontrary mean
ncp1=N*delta # non-central parameter
power.chi=pchisq(qchisq(1-alpha,df),df,ncp1,lower=FALSE)
q.alpha=qchisq(1-alpha,df) # 3.841
qcrit=qchisq(max(p.val,1-p.val),df,ncp=0)
# qcrit = prop.test(c(r1,r2),c(n1,n2))$statistic
curve(dchisq(x,df),0,2*qcrit,xlab=quote(chi),ylab='Density',bty='n',lwd=2) # H0
curve(dchisq(x,df,ncp1),add=TRUE,lty=2,col=2,lwd.2) # H1
# power.chi=AUC of right half of the H1 curve
abline(v=qcrit)
abline(v=q.alpha,col='gray',lty=3)
legend('topright',bty='n',legend=c('H0','H1','chisq@p.value','chisq@alpha'),col=c(1,2,1,'gray'),lty=c(1,2,1,3),lwd=c(2,2,1,1))
text(qcrit,0,round(qcrit,2))
y0=dchisq(qcrit,df,0)
y1=dchisq(qcrit,df,ncp=ncp1)
FPR.equal=y0/(y0+y1) # length ratio
FPR.less=p.val/(p.val+power.chi) # area ratio
FPR.alpha=alpha/(alpha+power.chi) # FPR before analysis
VAL=c(power=power.chi,p.value=p.val,FPR.equal=FPR.equal,
FPR.less=FPR.less,FPR.alpha=FPR.alpha)
print(VAL,digits=3)
invisible(VAL)
}
532:卵の名無しさん
18/05/19 07:58:23.94 rrC4yXIM.net
PPV2prevalence <- function(sensitivity,specificity,PPV) {
(1-specificity)*PPV/((1-specificity)*PPV)+(1-PPV)*sensitivity))
}
FPP2prior <- function(power,FPR=0.05,pval=0.05){
pval*(1-FPR)/(pval*(1-FPR)+FPR*power)
}
533:卵の名無しさん
18/05/19 10:40:37.24 oTDRH91u.net
PPV2prevalence <- function(sensitivity,specificity,PPV) {
(1-specificity)*PPV/((1-specificity)*PPV+(1-PPV)*sensitivity)
}
PPV2prevalence(0.75,0.99,0.9)
FPP2prior <- function(power,FPR=0.05,pval=0.05){
pval*(1-FPR)/(pval*(1-FPR)+FPR*power)
}
534:卵の名無しさん
18/05/20 11:45:44.17 n2fbjQMc.net
# These date were used in 1908 by W. S. Gosset ('Student')
# as an example to illustrate the use of his t test,
# in the paper in which the test was introduced.
A=c(0.7,-
535:1.6,-0.2,-1.2,-0.1,3.4,3.7,0.8,0.0,2.0) B=c(1.9,0.8,1.1,0.1,-0.1,4.4,5.5,1.6,4.6,3.4) t.test(A,B,var.equal = TRUE) mean(A) ; mean(B) sd(A) ; sd(B) (E=mean(B)-mean(A)) nA=length(A) ; nB=length(B) # The pooled estimate of the error within groups SEpooled=sqrt(weighted.mean(c(var(A),var(B)),c(nA-1,nB-1))) # Standard deviation of effect size SE=sqrt(SEpooled^2/nA+SEpooled^2/nB) (t.statistic=E/SE) 2*pt(t.statistic,nA+nB-2,lower=FALSE) # two-sided
536:卵の名無しさん
18/05/20 18:46:26.99 n2fbjQMc.net
# URLﾘﾝｸ(www.nejm.org)
r1=14
n1=840
r2=73 # placebo
n2=829
calc.FPR.p2(r1,r2,n1,n2)
calc.FPR0.p2(r1,r2,n1,n2)
prop.test(c(r1,r2),c(n1,n2))
prior.needed <- function(r1,r2,n1,n2,FPR=0.05){
pval=prop.test(c(r1,r2),c(n1,n2))$p.value
ES=pwr::ES.h(r1/n1,r2/n2)
power=pwr::pwr.2p2n.test(ES,n1,n2,sig.level=pval)$power
prior=pval*(1-FPR)/(pval*(1-FPR)+FPR*power)
return(prior)
}
prior.needed(r1,r2,n1,n2,FPR=0.05)
537:卵の名無しさん
18/05/21 00:13:06.30 r/Pmsrg8.net
# URLﾘﾝｸ(www.ncbi.nlm.nih.gov)
# p.508
# When p is a p-value with n1 samples, 95% ci of p-value of next experiment with n2 samples
# is supposed to be estimated by p2ci-function below.
p2ci <- function(p,n1=100,n2=100,sig.level=0.95){
lwr=pnorm(qnorm(p)*sqrt(n2/n1)-qnorm(1-(1-sig.level)/2)*sqrt(1+n2/n1))
# pnorm(qnorm(p1)-2.771808) # when n1=n2, 1.96 *s qrt(2) = 2.77
upr=pnorm(qnorm(p)*sqrt(n2/n1)+qnorm(1-(1-sig.level)/2)*sqrt(1+n2/n1))
# pnorm(qnorm(p1)+2.771808)
c(lwr,upr)
}
p2ci(0.05)
p2ci(0.001)
graphics.off()
pp=seq(0,1,length=1001)
plot(pp,sapply(pp,function(p)p2ci(p)[1]),type='l',bty='n',
xlab='initial p-value',ylab='C.I. of next p-value')
lines(pp,sapply(pp,function(p)p2ci(p)[2]))
###
(opt=optimise(function(p)p2ci(p)[1],c(0,1)))
p2ci(opt$minimum)
f <- function(n1=10,n2=10)c(t.test(rnorm(n1))$p.value,t.test(rnorm(n1))$p.value)
re=replicate(10^4,f())
points(re[1,],re[2,],col=rgb(0.01,0.01,0.01,0.01))
538:卵の名無しさん
18/05/22 15:20:40.05 pgTq+n72.net
等分散を仮定しないWelch法でのｔ検定からでも算出できるように
スクリプトを変更。
# URLﾘﾝｸ(www.physiology.org)
# Explorations in statistics: statistical facets of reproducibility Douglas Curran-Everett
p2p2w <- function(p.value,n1=10,n2=10,sd1=1,sd2=1,alpha=0.05){
p=p.value/2 # two.sided comparison
# Z-test
z=qnorm(p)
d.z=z*sqrt(sd1^2/n1+sd2^2/n2) # estimated difference of means by z.test
p2.z=pnorm(-qnorm(alpha/2)+z,lower=FALSE)
# Student
df=n1-1+n2-1
t=qt(p,df)
d.t=t*sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
p2.t=pt(-qt(alpha/2,df)+t,df,lower=FALSE)
# Welch
var1=sd1^2 ; var2=sd2^2
dfw=(var1/n1+var2/n2)^2 / (var1^2/n1^2/(n1-1)+var2^2/n2^2/(n2-1))
t.w=qt(p,dfw)
d.w=t.w*sqrt(var1/n1+var2/n2)
p2.w=pt(-qt(alpha/2,dfw)+t.w,dfw,lower=FALSE)
data.frame(p2.w, p2.t, p2.z)
}
539:卵の名無しさん
18/05/24 21:03:43.35 R1FtDNpg.net
# simulation on the assumption of the event ~ poisson distribution
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
stanString='
data{
int<lower=0> r1;
int<lower=0> r2;
int<lower=0> n1;
int<lower=0> n2;
}
parameters{
real<lower=0> lambda1;
real<lower=0> lambda2;
}
model{
r1 ~ poisson(lambda1) T[0,n1];
r2 ~ poisson(lambda2) T[0,n2];
}
'
# execute only on first occa
540:tion # stanModel=stan_model(model_code = stanString) # saveRDS(stanModel,'2x2p.rds')
541:卵の名無しさん
18/05/24 21:05:13.78 R1FtDNpg.net
prop.poisson <- function(r1,r2,n1,n2, alpha=0.05){
data <- list(r1=r1,r2=r2,n1=n1,n2=n2)
stanModel=readRDS('2x2p.rds')
fit=sampling(stanModel,data=data,chain=1,iter=10000)
print(fit)
ms=rstan::extract(fit)
k=length(ms$lambda1)
p.sim=numeric(k)
for(i in 1:k){
p.sim[i]=prop.test(c(ms$lambda1[i],ms$lambda2[i]),c(n1,n2))$p.value
}
BEST::plotPost(p.sim)
mean(p.sim<alpha)
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',prop.test(c(r1,r2),c(n1,n2))$p.value,'\n')
}
542:卵の名無しさん
18/05/25 22:02:04.22 8Ot6XASs.net
URLﾘﾝｸ(www.ncbi.nlm.nih.gov)
543:卵の名無しさん
18/05/30 10:36:17.99 5L7BTTj3.net
r=3
f<- function(x) (1-x)^(r-1)*x
curve (f(x))
auc=integrate (f,0,1)$value
pdf <- function (x) f(x)/auc
integrate(pdf, 0.5,1)$value
integrate (function (x)x*pdf(x),0,1)$value
2/(r+2)
z=3;N=9
f <- function (x) choose(N,z)*x^z*(1-x)^(N-z)
curve (f(x))
auc=integrate (f,0,1)$value
pdf <- function (x) f(x)/auc
integrate(pdf, 0.5,1)$value
integrate (function (x)x*pdf(x),0,1)$value
(z+1)/(N+2)
544:卵の名無しさん
18/05/30 20:54:56.95 w7KIivv+.net
f <- function(x) 1/sqrt(x*(1-x))
z <- function(x) integrate(f,0,x)$value
z(0.1)
z(0.5)
curve(f(x))
xx=seq(0,1,0.01)
plot(xx,sapply(xx,z))
545:卵の名無しさん
18/05/30 21:02:47.29 w7KIivv+.net
f <- function(x) 1/sqrt(x*(1-x))
z <- function(x) integrate(f,0,x)$value/pi
z(0.1)
z(0.5)
curve(f(x))
x=seq(.001,.999,by=.001)
plot(x,sapply(x,z))
abline(a=0,b=1)
546:卵の名無しさん
18/05/30 21:11:44.87 w7KIivv+.net
curve(pbeta(x,0.5,0.5))
547:卵の名無しさん
18/05/30 21:31:43.66 w7KIivv+.net
f <- function(x)dbeta(x,0.5,0.5)
zz <- function(x) integrate(f,0,x)$value
x=seq(.001,.999,by=.001)
plot(x,sapply(x,zz))
abline(a=0,b=1)
548:卵の名無しさん
18/06/01 11:46:08.93 UJHj2xQL.net
# marginal likelihoo of Fixed Effect Model
m.FE <- function(r1,n1,r2,n2,shape1=1,shape2=1){
choose(n1,r1)*choose(n2,r2)*beta(r1+r2+shape1,n1-r1+n2-r2+shape2)/beta(shape1,shape2)
}
# marginal likelihoo of Random Effect Model θ1~B(a1,ba),θ2~B(a2,b2)
m.REJ <- function(r1,n1,r2,n2,a1=1,b1=1,a2=1,b2=1){
2*choose(n1,r1)*choose(n2,r2)*integrate(
function(x){
dbeta(x,r1+a1,n1-r1+b1)*beta(r1+a1,n1-r1+b1)/beta(a1,b1)*
pbeta(x,r2+a2,n2-r2+b2,lower=FALSE)*beta(r2+a2,n2-r2+b2)/beta(a2,b2)
}
,0,1)$value
}
549:卵の名無しさん
18/06/02 14:59:28.21 iq0fCXui.net
Haldane <- function(x) 1/(x*(1-x))
curve(Haldane(x),0,1)
f <- function(x) integrate(Haldane,0.5,x)$value
x=seq(0.001,0.999,by=0.001)
plot(x,sapply(x,f),col=3)
curve(log(x/(1-x)),add=T)
550:卵の名無しさん
18/06/04 15:43:10.18 UVDFwwqx.net
エビデンス=周辺尤度= p(D |M1)=
∫
p(D | θ,M1)p(θ |M1)dθ
551:卵の名無しさん
18/06/04 15:43:29.56 UVDFwwqx.net
エビデンス=周辺尤度
= p(D |M1)=
∫
p(D | θ,M1)p(θ |M1)dθ
552:卵の名無しさん
18/06/04 18:46:43.68 UVDFwwqx.net
marginal likelihood of your model Mx is given by
p(D |Mx)= p(ξ1)p(D | ξ1)+p(ξ2)p(D | ξ2)+p(ξ3)p(D | ξ3)
=0.6×0.001+0.3×0.002+0.1×0.003
=0.0015.
The marginal likelihood is computed
553:卵の名無しさん
18/06/04 21:45:19.03 UVDFwwqx.net
stanStringJ='
data{
int<lower=0> N1;
int<lower=0> N2;
real Y1[N1];
real Y2[N2];
}
parameters{
real mu;
real<lower=0> sigma;
real alpha;
}
transformed parameters{
real mu1;
real mu2;
real delta;
real<lower=0> precision;
delta = alpha/sigma;
mu1 = mu + alpha/2;
mu2 = mu - alpha/2;
precision = 1/(sigma^2);
}
model{
Y1 ~ normal(mu1,sigma);
Y2 ~ normal(mu2,sigma);
mu ~ uniform(-9999,9999);
precision ~ gamma(0.5,0.5);
delta ~ cauchy(0,1);
}
'
554:卵の名無しさん
18/06/04 21:46:02.48 UVDFwwqx.net
N1=n1=137 ; m1=23.8 ; sd1=10.4 ; N2=n2=278 ; m2=25.2 ; sd2=12.3
U=scale(rnorm(n1))*sd1+m1 ; E=scale(rnorm(n2))*sd2+m2
Y1=as.vector(U) ; Y2=as.vector(E)
data=list(N1=N1,N2=N2,Y1=Y1,Y2=Y2)
# execute only on first occation
# JZS.model=stan_model(model_code = stanStringJ)
# JZS.model=stan_model('JZS.stan')
# saveRDS(JZS.model
555:,'JZS.rds') JZS.model=readRDS('JZS.rds') fitJ=sampling(JZS.model,data=data, iter=20000) print(fitJ,probs = c(.025,.5,.975)) ms=rstan::extract(fitJ) dES=density(ms$delta) plot(dES,lwd=2, xlab='Effect Size (δ)',main='The Savage?Dickey method',bty='l') curve(dcauchy(x),lty=3, add=TRUE) abline(v=0,col=8) (d1=dES$y[which.min(dES$x^2)]) # posterior density at ES=0 (d0=dcauchy(0)) # prior density at ES=0 points(0,d0,cex=2) points(0,d1,pch=19,cex=2) legend('topleft',bty='n',legend=c('prior','posterior'),lty=c(3,1),lwd=c(1,2)) (BF10=d1/d0) text(0,0,paste('BF10=',round(BF10,2))) library(BayesFactor) 1/exp(ttestBF(Y1,Y2,rscale = 1)@bayesFactor[['bf']])
556:卵の名無しさん
18/06/05 09:10:30.95 Ob/ggqKh.net
N=10000
x=rbeta(N,0.5,0.5)
y=rbeta(N,0.5,0.5)
z=x-y
hist(z,freq=F)
dz=density(z)
dz$y[which.min(dz$x^2)]
z=x/y
hist(log10(z),freq=F)
dz=density(z)
min=1
dz$y[which.min((dz$x-min)^2)]
557:卵の名無しさん
18/06/05 09:14:46.94 Ob/ggqKh.net
data{ //binomBF.stan
int r1;
int r2;
int n1;
int n2;
real shape1;
real shape2;
}
parameters{
real <lower=0, upper=1> theta1;
real <lower=0, upper=1> theta2;
real <lower=0, upper=1> th1;
real <lower=0, upper=1> th2;
}
transformed parameters{
real rd;
real rr;
real rd0;
real rr0;
rd = theta1 - theta2;
rd0 = th1 - th2;
rr = theta1/theta2;
rr0 = th1/th2
}
model{
r1 ~ binomial(n1,theta1);
r2 ~ binomial(n2,theta2);
theta1 ~ beta(shape1,shape2);
theta2 ~ beta(shape1,shape2);
th1 ~ beta(shape1,shape2);
th2 ~ beta(shape1,shape2);
}
558:卵の名無しさん
18/06/06 21:19:25.17 IpiYYmjt.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)
559:卵の名無しさん
18/06/06 21:19:31.46 IpiYYmjt.net
一方、十八という数字が好きな太郎は一人ずつ調べて１８人めの女子がみつかったところで調査を終えることにした。
１８人めがみつかったのは花子と同じく５０人めであった。
帰無仮説のもとで
１８人がみつかるのが５０人めである確率は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)
主観である、検定の中止の基準の差でｐ値や信頼区間が変化するのは変だという批判である。
560:卵の名無しさん
18/06/07 07:36:58.84 6l5aJg03.net
50C18 * 0.5^18 * 0.5^32 �
561:ﾆ 49C17 * 0.5^17 * 0.5^32 * 0.5 の違いでしょう 18人目を見つけた人数を調べるというデザインがおかしいよね これ事前確率0.5で50人調査して女が18人っていうのを ベイズ更新していったらどうなる？
562:卵の名無しさん
18/06/07 07:47:39.40 Ni5wt/sw.net
>>524
酷いのになるとp<0.05になったらやめるとかいうのもあるな。
p-hackingと呼ばれる
563:卵の名無しさん
18/06/08 22:06:59.94 dTNUKiNw.net
r=21
N=20
a=0.5
b=0.5
p=a/(a+b+r+(1:N))
q=cumsum(p)
q
plot(1:N,p,ann=F)
plot(1:N,q,ann=F)
​
564:卵の名無しさん
18/06/12 00:02:23.51 XLL1LdWn.net
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=((z*1.0)/(N*1.0)-FP)/(TP-FP);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2);
}
')
565:卵の名無しさん
18/06/12 00:34:11.64 XLL1LdWn.net
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=((z*1.0)/(N*1.0)-FP)/(TP-FP);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2) T[0.5,];
}
')
model=stan_model(model_code = stanString)
fit=sampling(model,data=data,control=list(adapt_delta=0.99),iter=10000)
print(fit,digits=3,probs=c(.025,.50,.975))
566:卵の名無しさん
18/06/12 19:09:54.64 XLL1LdWn.net
# pdfからcdfの逆関数を作ってHDIを表示させて逆関数を返す
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
print(c(hdi[1],hdi[2]),digits=5)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}
567:卵の名無しさん
18/06/12 19:12:46.61 XLL1LdWn.net
.n=50
.r=40
b=0.01
f = function(x,a,b,n,r){ # x:prevalence, a:TP, b:FP=0.01
p=x*a+(1-x)*b
choose(n,r)*p^r*(1-p)^(n-r)
}
f2 = function(x,n=.n,r=.r){
cubature::adaptIntegrate(function(ab)f(x,ab[1],b,n,r),
c(0,0),c(1,1))$integral
}
vf2=Vectorize(function(x)f2(x,n=.n,r=.r))
curve(vf2(x),ylab='',yaxs='i',axes=FALSE,lwd=2,xlab='prevalence') ; axis(1)
# points(0:10/10,vf2(0:10/10),pch=19)
optimise(function(x) vf2(x),c(0,1),maximum = TRUE)
auc=integrate(function(x)vf2(x),0,1)$value
pdf<-function(x)vf2(x)/auc
vpdf=Vectorize(pdf)
integrate(function(x)x*pdf(x),0,1)$value # mean
cdf=function(x)integrate(pdf,0,x)$value
vcdf=Vectorize(cdf)
# time consuming processes
# curve(vcdf(x),bty='l',n=64)
inv_cdf <- function(u){
uniroot(function(x,u0=u)vcdf(x)-u0,c(0,1),tol = 1e-18)$root
}
vinv_cdf=Vectorize(inv_cdf)
# curve(vinv_cdf(x),lwd=2,n=64)
hdi=HDInterval::hdi(vinv_cdf) ; hdi
# lower upper
# 0.7456652 1.0000000
568:卵の名無しさん
18/06/12 19:23:42.87 XLL1LdWn.net
# random numbers following PDF by John von Neuman's method
vonNeumann2 <- function(PDF,xmin
569:=0,xmax=1,N=10000,Print=TRUE,...){ xx=seq(xmin,xmax,length=N+1) xx=xx[-(N+1)] xx=xx[-1] ymax=max(PDF(xx)) Ux=runif(N,xmin,xmax) Uy=runif(N,0,ymax) Rand=Ux[which(Uy<=PDF(Ux))] if(Print){ hist(Rand,xlim=c(xmin,xmax),freq=FALSE,col=sample(colors(),1),main='',...) AUC=integrate(PDF,xmin,xmax)$value lines(xx,sapply(xx,function(x)PDF(x)/AUC)) } hdi=HDInterval::hdi(Rand) print(c(hdi[1],hdi[2]),digits=5) invisible(Rand) }
570:卵の名無しさん
18/06/12 20:53:05.85 Ex3k8fq/.net
>>531
ここでもうりゅう先輩が迷惑掛けてんのか？
ウリュウなあ
こいつはなあ、生まれついてのビッグマウスであちこちに自分を売り込むが、
卒業しても国試浪人で医師免許ない50過ぎでは相手にされない
国試対策塾で非常識講師で細々と食つなぐが学生に馬鹿にされる
自分の医師コンプを隠すために医学生たちを「底辺」などという
実は自分が凄まじい底辺なのだが気づいていない
こんな嘘つきデブがのさばっているスレだな
ご苦労なこったよ、うりゅうのおっさん
わからねえとでも思ってんだろどーせ
571:卵の名無しさん
18/06/13 08:53:38.73 rcG4xYvl.net
library(rjags)
N=50
z=40
FP=0.01
shape1=1
shape2=1
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("TP","x","theta"), n.iter=100000 )
js=as.matrix(codaSamples)
head(js)
BEST::plotPost(js[,'TP'],xlab='sensitivity')
BEST::plotPost(js[,'x'],xlab='prevalence')
BEST::plotPost(js[,'theta'],xlab='positive result',showMode = TRUE)
572:卵の名無しさん
18/06/16 16:23:17.07 V7A3qxKV.net
経理課の須藤は着服をやめろ！
勤務実態もないのに、グループ病院内から管理手当て（10万円）をもらうな！！！
意図的な給与操作、どうにかしろ！
573:卵の名無しさん
18/06/17 05:51:35.25 ifh2AARM.net
seqN <- function(N=100,K=5){
a=numeric(N)
for(i in 1:K) a[i]=2^(i-1)
for(i in K:(N-1)){
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=a[i+1]+a[i-j] # recursion formula
}
}
P0=numeric(N)
for(i in 1:N) P0[i]=a[i]/2^i # P0(n)=a(n)/2^n
P0
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,1]=P0
head(MP);tail(MP)
MP[1,2]=1/2
for(i in (K-2):K) MP[1,i]=0
for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=1/2*MP[i,k-1]
} # Pk(n+1)=1/2*P(k-1)(n)
ret=1-apply(MP,1,sum)
ret[N]
}
seqN(100,5)
seqN(1000,10)
574:卵の名無しさん
18/06/17 07:58:05.56 ifh2AARM.net
## 表の出る確率がｐであるとき、N回コインを投げて K回以上表が連続する確率
seqNp <- function(N=100,K=5,p=0.5){
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(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
575: } P0=numeric(N) for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n P0 MP=matrix(rep(NA,N*K),ncol=K) colnames(MP)=paste0('P',0:(K-1)) MP[,'P0']=P0 head(MP);tail(MP) MP[1,'P1']=p for(i in (K-2):K) MP[1,i]=0 for(k in 2:K){ for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1] } # Pk(n+1)=p*P(k-1)(n) ret=1-apply(MP,1,sum) ret[N] }
576:卵の名無しさん
18/06/17 08:33:04.03 LIEnVEKd.net
p:表
q=1-p
# Pk(n) (k=0,1,2,3,4)を途中、5連続して表が出ていなくて
# 最後のk回は連続して表が出ている確率とする。
#
P0(1)=q
P1(1)=p
P2(1)=P3(1)=P4(1)=0
P(k+1)(n+1)=p*Pk(n)
P0(n+1)=q*{P0(n)+P1(n)+P2(n)+P3(n)+P4(n)}
=q*{P0(n)+p*P0(n-1)+p^2*P0(n-2)+p^3*P0(n-3)+p^4*P0(n-4)}
P0(n)=a(n)*p^n
# a(n+1)p^(n+1)=q*p^n{a(n)+a(n-1)+a(n-2)+a(n-3)+a(n-4)}
# a(n+1)=q/p1*(a(n)+a(n-1)+a(n-2)+a(n-3)+a(n-4))
a(n)=P0(n)/p^n
577:卵の名無しさん
18/06/17 09:43:49.13 ifh2AARM.net
>>532
統計くらいできるのが国立卒の普通の臨床医。
おい、ド底辺
統計処理からはおまえは
都外のド底辺シリツ医大卒と推測されたが、あってるか？
578:卵の名無しさん
18/06/17 09:45:16.20 ifh2AARM.net
## 表の出る確率がｐであるとき、N回コインを投げて K回以上表が連続する確率に一般化してみた。
seqNp <- function(N=100,K=5,p=0.5){
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(i in K:(N-1)){ # recursive formula
a[i+1]=0
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)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
head(MP);tail(MP)
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0
for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)
ret[N]
}
579:卵の名無しさん
18/06/17 22:30:10.45 ifh2AARM.net
# pdfからcdfの逆関数を作ってHDIを表示させて逆関数を返す
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
print(c(hdi[1],hdi[2]),digits=5)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}
580:卵の名無しさん
18/06/18 20:41:20.64 G5tabFnb.net
library(rjags)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
N=100 ; FP=0.01 ; shape1=1 ; shape2=1
guess.TP <- function(z){
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples(jagsModel,variable=c("TP","x"), n.iter=10000)
js=as.matrix(codaSamples)
m.TP=mean(js[,'TP'])
ci.TP=HPDinterval(as.mcmc(js[,'TP']))
m.x=mean(js[,'x'])
ci.x=HPDinterval(as.mcmc(js[,'x']))
c(m.TP=m.TP,ci.TP=ci.TP,m.x=m.x,ci.x=ci.x)
}
zz=1:20*5
re=sapply(zz,guess.TP)
head(re[,1:4])
re=as.matrix(re)
plot(zz,re['m.TP',],bty='l',ylim=c(0,1),type='n',las=1,
xlab='n : positives out of 100',ylab='sensitivity')
segments(zz,re[2,],zz,re[3,],col=8,lwd=3)
points(zz,re['m.TP',],pch=16)
581:卵の名無しさん
18/06/18 22:45:56.48 G5tabFnb.net
guess.TP2 <- function(z,FP){
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples(jagsModel,variable=c("TP"), n.iter=10000,thin=5)
js=as.matrix(codaSamples)
mean(js[,'TP'])
}
vG=Vectorize(guess.TP2)
n=1:20*5
FP=seq(0,0.25,length=20)
TP=outer(n,FP,vG) # wait several minutes
contour(n,FP,TP, col='navy',
xlab='n : positives out of 100',ylab='FP : 1-specificity',bty='l',nlevels=64)
points(50,0.10,pch='+',col='red',cex=1.5)
582:卵の名無しさん
18/06/19 22:41:00.16 ant1u+bV.net
n=5 # prime number
nn=1:(n-1)
tasu <- function(x,y) (x+y)%%n
hiku <- function(x,y) (x-y)%%n # row - col
kake <- function(x,y) (x*y)%%n
g=function(x) nn[which(x==1)]
.M=outer(nn,nn,kake)
G=apply(.M,2,g)
gyaku <- function(x) nn[which(G==(x%%n))]
waru <- function(x,y) (x*gyaku(y))%%n # row / col
waru(3,2)
xx=yy=c(0,nn)
names(xx)=paste0('x',c(0,nn))
names(yy)=paste0('y',c(0,nn))
outer(xx,yy,tasu) # x + y
outer(xx,yy,hiku) # x - y
outer(xx,yy,kake) # x * y
X=Y=nn
outer(X,Y,waru) # WRONG!!
outer(X,Y,Vectorize(waru))
a=expand.grid(X,Y)
b=matrix(mapply(waru,a[,1],a[,2]),ncol=length(X))
rownames(b)=paste0('x',nn)
colnames(b)=paste0('y',nn)
b # x / y
583:卵の名無しさん
18/06/20 19:25:59.38 myhyWcyK.net
rule3 <- function(n,confidence.level=0.95){
p=1/n
q=1-p # q^n.sample > 1-confidence.level
n.sample = log(1-confidence.level)/log(q)
return(n.sample)
}
n.sample=log(0.05)/log(0.999)
584:卵の名無しさん
18/06/22 14:47:07.52 ETsHYxXe.net
shuffle <- function(Cards){
n=length(Cards)
n1=ceiling(n/2)
n2=n-n1
C1=Cards[1:n1]
C2=Cards[(n1+1):n]
ret=NULL
for(i in 1:n1){
ret=c(ret,C1[i],C2[i])
}
ret[is.na(ret)==F]
}
x=as.character(c('A',2:10,'J','Q','K'))
cat(x,'\n') ; cat(shuffle(x))
Shuffles <- function(x){
tmp=shuffle(x)
i=1
while(!identical(x,tmp)){
tmp=shuffle(tmp)
i=i+1
}
return(i)
}
f =function(x)Shuffles(1:x)
nn=1:53
y=sapply(nn,f)
plot(nn,y,pch=16,bty='l',xlab='cards',ylab='shuffles')
cbind(nn,y)
585:卵の名無しさん
18/06/22 20:07:08.88 ETsHYxXe.net
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)
}
x=c(4, 3, 5, 2, 1)
inversion(x)
is.even= function(x) !inversion(x)%%2
is.even(x)
prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx]) ; is.even(X)
Y=numeric(n-1)
for (i in 1:(n-1)){
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # 囚人iが見えない番号
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()))
586:卵の名無しさん
18/06/26 07:03:32.73 kT/81/Gu.net
# 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()))
587:卵の名無しさん
18/06/26 10:42:11.82 kT/81/Gu.net
inversion <- function(x){ #inversion number
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret)
}
is.even <- function(x) !inversion(x)%%2 # even inversion?
even.perm <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx]) # row of 99 prisoner numbers
is.even(X)
}
mean(replicate(1e3,even.perm())) # probability of even permutation
588:卵の名無しさん
18/06/27 14:01:31.44 gge/PUDl.net
#ある大学の学生数は５００以上１０００人以下であることはわかっている。
#無作為に２人を抽出して調べたところ
#二人とも女子学生である確率は１／２であった。
#この大学の女子学生数と男子学生数は何人か？
girlsboys <- function(g,b) g*(g-1)/(g+b)/(g+b-1)==1/2
gr=expand.grid(1:1000,1:1000)
(re=gr[which(mapply(girlsboys,gr[,1],gr[,2])),])
girlsboys(re[nrow(re),1],re[nrow(re),2])
589:卵の名無しさん
18/06/27 14:03:56.49 gge/PUDl.net
# ある大学の学生数は５００以上１０００人以下であることはわかっている。
# 無作為に２人を抽出して調べたところ
# 二人とも女子学生である確率は１／２であった。
# この大学の女子学生数と男子学生数は何人か？
girlsboys <- function(g,b) g*(g-1)/(g+b)/(g+b-1)==1/2
gr=expand.grid(1:1000,1:1000)
(re=gr[which(mapply(girlsboys,gr[,1],gr[,2])),])
# 検証
Vectorize(girlsboys)(re[,1],re[,2])
590:卵の名無しさん
18/06/27 17:53:09.08 gge/PUDl.net
N=120
r=10
D=c(rep(1,r),rep(0,N-r))
hiseiki <- function(m){
found=0
for(i in 1:(N-r+m)){
found=found+sample(D,1)
if(found==m) break
}
return(i)
}
re=replicate(1e4,hiseiki(3))
mean(re)
sd(re)
BEST::plotPost(re,breaks=30)
591:卵の名無しさん
18/06/29 07:33:02.38 zpGshx3p.net
cereals <- function(n=5){
coupons=NULL
while(!all((1:n) %in% coupons)){
coupons=append(sample(1:n,1),coupons)
}
return(length(coupons))
}
re=replicate(100,mean(replicate(1e3,cereals(5))))
> mean(re)
[1] 11.43503
592:卵の名無しさん
18/06/29 08:13:57.68 zpGshx3p.net
p=4:1
cereals <- function(p){
n=length(p)
coupons=NULL
while(!all((1:n) %in% coupons)){
coupons=append(sample(1:n,1,p=p),coupons)
}
return(length(coupons))
}
mean(replicate(1e3,cereals(p)))
re=replicate(100,mean(replicate(1e3,cereals(p))))
mean(re)
593:卵の名無しさん
18/06/29 16:00:32.65 zpGshx3p.net
blood.type <- function(p,need){
n=length(p)
ABO=NULL
enough <- function(x){
pool=numeric(n)
for(i in 1:n) pool[i]=sum(ABO==i)
all(pool >= need)
}
while(!enough(ABO)){
ABO=append(sample(1:n,1,p=p),ABO)
}
return(length(ABO))
}
p=4:1
need=c(10,10,5,2)
re=replicate(1e4,blood.type(p,need))
BEST::plotPost(re)
594:卵の名無しさん
18/06/30 05:47:26.70 BDmYWstD.net
BT <- function (a,b,c,d) 1/a + 1/b + 1/c + 1/d - 1/(a+b) - 1/(a+c) - 1/(b+c) - 1/(a+d) - 1/(b+d) - 1/(c+d) + 1/(a+b
595:+c) + 1/(d+a+b) + 1/(c+d+a) + 1/(b+c+d) - 1/(a+b+c+d) a=4 b=3 c=2 d=1 s =a+b+c+d a=a/s b=b/s c=c/s d=d/s BT(a,b,c,d)
596:卵の名無しさん
18/06/30 13:48:31.26 j2CU1Lw0.net
even.tally <- function(a=3 , b=2){
idx=combn(1:(a+b),a)
n=ncol(idx)
mat=matrix(0,nrow=n,ncol=a+b)
for(i in 1:n) mat[i,idx[,i]]=1
tally <- function (x) any(cumsum(x)==cumsum(1-x))
mean(apply (mat,1,tally))
}
even.tally()
even.tally(5,10)
597:卵の名無しさん
18/06/30 15:09:55.32 j2CU1Lw0.net
a=750 ; b=250
v=c(rep(1,a),rep(0,b))
f <- function(v){
x=sample(v)
any(cumsum(x)==cumsum(1-x))
}
mean(replicate(1e5,f(v)))
598:卵の名無しさん
18/07/01 10:28:46.39 F3KIhVUK.net
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x)
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}
#
birth(100, 3)
vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
599:卵の名無しさん
18/07/01 11:53:44.57 F3KIhVUK.net
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x))
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}
#
birth(100, 3)
vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
600:卵の名無しさん
18/07/01 12:03:30.64 F3KIhVUK.net
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x)
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}
#
birth(100, 3)
vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
601:卵の名無しさん
18/07/01 12:08:34.56 F3KIhVUK.net
インフルエンザの迅速キットは特異度は高いが感度は検査時期によって左右される。
ある診断キットが開発されたとする。
このキットは特異度は９９％と良好であったが、
感度については確かな情報がない。
事前確率分布として一様分布を仮定する。
５０人を無作為抽出してこの診断キットで診断したところ４０人が陽性であった。
この母集団の有病率の期待値と９５％信用区間はいくらか？
またこの診断キットの感度の期待値と９５％信用区間はいくらか
暇つぶしにこれをMCMCを使わずに解く方法を考えていた。
偽陽性率FP=0.01として
陽性確率p=TP*x+(1-x)*FP
尤度が50C40*p^40*(1-p)^10
TPは一様分布なので積分消去して
確率密度関数に比例する関数を作ってarea under the curveで割って確率密度関数化したのち積分して累積密度関数をつくる。この累積密度関数の逆関数を作って95%区間が最短になる区間を計算すれば信頼区間が算出できる。
この結果がstanでのシミュレーションの結果と一致すればよし。
602:卵の名無しさん
18/07/01 22:00:46.73 lSk1MYzX.net
# # choose(n,r) == gamma(n+1) / (gamma(r+1) * gamma(n-r+1))
same.birthday <- function(n) 1-choose(365+97/400,n)*factorial(n)/(365+97/400)^n
plot(x,y,bty='l',xlab='subjects',ylab='probability')
curve(same.birthday(x),add = TRUE)
abline(h=0.5,col=8)
same.birthday(22:23)
603:卵の名無しさん
18/07/02 01:26:48.74 pp47QgIN.net
library(rjags)
N=50
z=40
FP=0.01
shape1=1
shape2=1
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("TP","x","theta"), n.iter=100000 )
js=as.matrix(codaSamples)
head(js)
BEST::plotPost(js[,'TP'],xlab='sensitivity')
BEST::plotPost(js[,'x'],xlab='prevalence')
BEST::plotPost(js[,'theta'],xlab='positive result',showMode = TRUE)
604:卵の名無しさん
18/07/02 01:28:46.92 pp47QgIN.net
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
real<lower=0,upper=1> x; //prevalence
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=TP*x+FP*(1-x);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2); // T[0.5,];
x ~ beta(shape1,shape2);
}
')
605:卵の名無しさん
18/07/02 01:29:10.57 pp47QgIN.net
# model=stan_model(model_code = stanString)
# saveRDS(model,'quick_kit.rds')
model=readRDS('quick_kit.rds')
fit=sampling(model,data=data,iter=10000)
print(fit,digits=3,probs=c(.025,.50,.975))
stan_trace(fit)
# stan_diag(fit)
stan_ac(fit)
stan_dens(fit,separate_chains = TRUE)
stan_hist(fit,fill='skyblue',bins=15,pars=c('x','TP'))
ms=rstan::extract(fit)
BEST::plotPost(ms$TP,showMode = TRUE,xlab='TP')
BEST::plotPost(ms$x,showMode = FALSE,xlab='prevalence',col=sample(colours(),1))
606:卵の名無しさん
18/07/02 21:40:12.88 tMXA02ZR.net
escape.cliff <- function(p=2/3,k=1000){
pos=1
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1)
607:,1,prob=c(p,1-p)) } return(TRUE) } mean(replicate(1e2,escape.cliff()))
608:卵の名無しさん
18/07/02 22:58:14.27 tMXA02ZR.net
escape.cliff <- function(pos=1,p=2/3,k=1000){
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e3,escape.cliff(10,0.5,100)))
#
totter <- function(p=2/3,pos=1,pos0=NULL,k=1e3){
if(is.null(pos0)) pos0=pos
for(i in 1:k) {
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(pos>pos0)
}
mean(replicate(1e3,totter(0.5,10,5 )))
609:卵の名無しさん
18/07/03 07:46:06.00 22JZXDLY.net
exam <- function(p=0.5,hit=0,money=5,k=5){
while(hit<k){
money=money-1
if(money==0) return(FALSE)
shoot=sample(c(1,0),1,prob=c(p,1-p))
if(shoot){
money=money+1
hit=hit +1
}
}
return(TRUE)
}
mean(replicate(1e5, exam()))
610:卵の名無しさん
18/07/03 08:47:34.93 ltuSOSv2.net
ryunen <- function(p=0.6 ,money=10,grade=0,ryu=0){
while(grade<6){
test=rbinom(1,1,p)
if(test) grade=grade+1
else ryu=ryu+1
if(ryu==2){
money=money-5
ryu=ryu-1
}
if(money<=0) return(FALSE)
}
return(grade==6)
}
mean(replicate(1e5,ryunen()))
611:卵の名無しさん
18/07/03 11:46:57.62 22JZXDLY.net
escape.cliff <- function(pos=1,p=2/3,k=1000){
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e5,escape.cliff(1,2/3,3)))
4/7
612:卵の名無しさん
18/07/03 13:41:51.20 22JZXDLY.net
gmbl <- function(money=20,p=18/38){
while(0<money & money <40){
money=money+sample(c(1,-1),1,p=c(p,1-p))
}
return(money==40)
}
mean(replicate(1e3, gmbl()))
613:卵の名無しさん
18/07/03 14:41:00.85 22JZXDLY.net
x=c(rep(1,4),rep(0,48))
mean(replicate(1e5,which(sample(x)==1)[1]))
y=numeric(48)
for(i in 1:48) y[i]=i*choose(48,i-1)*factorial(i-1)*4*factorial(52-i)/factorial(52)
sum(y)
614:卵の名無しさん
18/07/03 14:42:24.05 22JZXDLY.net
(a) A railroad numbers its locomotives in order, 1, 2, . . . , N. One day you see a locomotive and its number is 60.
Guess how many locomotives the company has.
(b) You have looked at 5 locomotives and the largest number observed is 60.
Again guess how many locomotives the company has.
615:卵の名無しさん
18/07/03 16:06:53.54 22JZXDLY.net
n=5
m=60
N=n:100
loco <- function(x){
max(sample(1:x,n))
}
vloco=Vectorize(loco)
loco.sim <- function() N[which(vloco(N)==m)]
locomotives=unlist(replicate(1e4,loco.sim()))
BEST::plotPost(locomotives,breaks=30)
HDInterval::hdi(locomotives)
616:卵の名無しさん
18/07/03 18:26:54.88 22JZXDLY.net
n=5
library(gtools)
y=permutations(n,n,1:n)
f=function(x) sum(x==1:n)
sum(apply(y,1,f))/length(y)
1/n
617:卵の名無しさん
18/07/03 18:39:40.14 22JZXDLY.net
f = function(n){
sum(sample(1:n)==1:n)/n
}
g=function(n)mean(replicate(1e4,f(n)))
vg=Vectorize(g)
x=1:100
plot(x,vg(x))
curve(1/x,add=T)
618:卵の名無しさん
18/07/03 21:26:39.83 22JZXDLY.net
gr=expand.grid(1:9,1:9,1:9)
u=10^(2:0)
f=function(a,b,c,x=1776) sum(c(a,b,c)*u+c(b,c,a)*u+c(c,a,b)*u)-x
# when A>B>C
idx=which(mapply(f,gr[,1],gr[,2],gr[,3])==0 & gr[,1]>gr[,2] & gr[,2]>gr[,3])
gr[idx,]
# otherwise
idx2=which(mapply(f,gr[,1],gr[,2],gr[,3])==0)
length (idx2)
gr[idx2,]
619:卵の名無しさん
18/07/04 00:37:12.00 sDwIwpA3.net
n=5
m=60
N=n:100
loco <- function(x){
max(sample(1:x,n))
}
vloco=Vectorize(loco)
loco.sim <- function() N[which(vloco(N)==m)]
locomotives=unlist(replicate(1e4,loco.sim()))
BEST::plotPost(locomotives,breaks=30)
HDInterval::hdi(locomotives)
n=60:100
pmf=choose(59,4)/choose(n,5) #Pr(max=60|n)
pdf=pmf/sum (pmf)
sum( n*pdf) #E(n)
lines(n,pdf)
620:卵の名無しさん
18/07/05 10:51:42.59 9YHvFei/.net
own <- function(n){
mean(replicate(1e5,sum(sample(1:n)-1:n==0)))
}
own(100)
621:卵の名無しさん
18/07/05 10:54:03.43 9YHvFei/.net
r1 <- function(x){ # rotate by one bead
n=length(x)
y=numeric(n)
y[1]=x[n]
for(i in 1:(n-1)){
y[i+1]=x[i]
}
return(y)
}
rn <- function(x){ # every rotation
n=length(x)
mat=matrix(rep(NA,n^2),ncol=n)
mat[,1]=x
for(i in 1:(n-1)){
mat[,i+1]=r1(mat[,i])
}
return(t(mat))
}
same <- function(x,y){
if(sum(x)!=sum(y)) return(FALSE)
f=function(a,b=y){ # is equal to y
all(a==b)
}
mat=rbind(rn(x),rn(rev(x))) # with symmetric conversion
any(apply(mat,1,f))
}
622:卵の名無しさん
18/07/05 10:55:44.62 9YHvFei/.net
dec2bin <- function(num, digit=0){ # decimal to 0,1 vector
if(num <= 0 && digit <= 0){
return(NULL)
}else{
return(
623:append(Recall(num%/%2,digit-1), num%%2)) } } vd2b=Vectorize(dec2bin) # bracelett <- function(n){ mat=t(vd2b(0:(2^n-1),n)) # make all permutation of beads # head(mat) ret=list() # list of the same bracelett for(i in 1:2^n){ ret[[i]]=which(apply(mat,1,function(z)same(z,mat[i,]))) } # head(ret) ; table(unlist(ret)) del=NULL for(i in 1:2^n) del=append(del,ret[[i]][-1]) 2^n - length(unique(del)) }
624:卵の名無しさん
18/07/05 14:59:01.13 bn3tqJqU.net
> library(gtools)
> swap <- function(n){
+ perm=permutations(n,n,1:n)
+ mean(apply(perm,1,function(x)sum(x==1:n)))
+ }
> swap(9)
[1] 1
> swap.sim <- function(n,k=1e5){
+ mean(replicate(k,sum(sample(1:n)==1:n)))
+ }
> swap.sim(9)
[1] 1.0004
625:卵の名無しさん
18/07/05 20:18:51.37 bn3tqJqU.net
swap2 <- function(n){ # k*nCk*(1/n)^k*(1-1/n)^(n-k) k=0:n
re=numeric(n)
for(i in 0:n){
p=1/n
re[i]=i*choose(n,i)*p^i*(1-p)^(n-i)
}
sum(re)
}
swap2(7)
#
swap3 <- function(n){
re=numeric(n)
for(i in 0:n){
p=1/n
re[i]=i*dbinom(i,n,p)
}
sum(re)
}
swap3(7)
626:卵の名無しさん
18/07/06 09:11:30.20 UEcW6fma.net
f012 <- function(n){
args=list()
length(args)=n-1
args[[1]]=args[[n-1]]=1:2
for(i in 2:(n-2)){
args[[i]]=0:2
}
gr=do.call(expand.grid,args)
gr=as.matrix(gr)
ret=numeric()
for(i in 1:nrow(gr)){
ret[i]=all(diff(gr[i,])!=0)
}
sum(ret)
}
627:卵の名無しさん
18/07/21 16:31:14.15 V1Wm3iKf.net
/* SEND+MORE=MONEYの覆面暗算をC言語で解いてみる。*/
#include<stdio.h> /* SEND+MORE=MONEY */
int compare_int(const void *a, const void *b)
{return *(int*)a - *(int*)b;}
int i,j;
int unique(int num[8]){
qsort(num,8,sizeof(int),compare_int);
for(i=0;i<8;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;
}
main(){
int S,E,N,D,M,O,R,Y;
for(S = 1; S < 10; S++){
for(E = 0; E < 10; E++){
for(N = 0; N < 10; N++){
for(D = 0; D < 10; D++){
for(M = 1; M < 10; M++){
for(O = 0; O < 10; O++){
for(R = 0; R < 10; R++){
for(Y = 0; Y < 10; Y++){
if(S*1000+E*100+N*10+D+M*1000+O*100+R*10+E==M*10000+O*1000+N*100+E*10+Y){
int num[]={S,E,N,D,M,O,R,Y};
if(unique(num)==1){
printf("%d%d%d%d+%d%d%d%d=%d%d%d%d%d\n", S,E,N,D,M,O,R,E,M,O,N,E,Y);
}}}}}}}}}}}
C:\MinGW>gcc sendmore.c -o money
C:\MinGW>money
9567+1085=10652
628:卵の名無しさん
18/07/22 11:46:07.83 tW9s1ZUi.net
#include<stdio.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;}
int unique(int num[]){
int i,j,n=10;
qsort(num,n,sizeof(int),compare_int);
for(i=0;i<n;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;}
main(){ /* ド底辺＋私立医＝裏口馬鹿 */
int n=1,A,B,C,D,E,F,G,H,I,J;
for(A = 1; A < 10; A++){
for(B = 0; B < 10; B++){
for(C = 0; C < 10; C++){
for(D = 1; D < 10; D++){
for(E = 0; E < 10; E++){
for(F = 0; F < 10; F++){
for(G = 1; G < 10; G++){
for(H = 0; H < 10; H++){
for(I = 0; I < 10; I++){
for(J = 0; J < 10; J++){
if(A*100+B*10+C +D*100+E*10+F==G*1000+H*100+I*10+J){
int num[]={A,B,C,D,E,F,G,H,I,J};
if(unique(num)==1){
printf("%2d: %d%d%d + %d%d%d = %d%d%d%d\n", n,A,B,C,D,E,F,G,H,I,J);
n++;
}}}}} }}} }}} }}
629:卵の名無しさん
18/07/23 16:11:37.39 B1Q5jsGg.net
# ２つの整数があります。
# それらをたしてできた数は、
# 十の位と一の位の数字が等しい２けたの整数になり、
# それらをかけてできた数は、
# 百の位、十の位、一の位が等しい３けたの整数になりました。
# このような２つの整数の組をあるだけ答えなさい。
gr=expand.grid(1:99,1:99)
f <- function(x,y){ all(x>=y,c((x+y)%%10==(x+y)%/%10,(x*y)%/%100==((x*y)%/%10)%%10,(x*y)%/%100==(x*y)%%10))}
gr[which(mapply(f,gr[,1],gr[,2])==TRUE),]
630:卵の名無しさん
18/07/23 17:34:03.68 B1Q5jsGg.net
# 実行速度遅すぎて実用性なし
# Consider a function which, for a given whole number n,
# returns the number of ones required when writing out all numbers between 0 and n.
# For example, f(13)=6. Notice that f(1)=1.
# What is the next largest number n such that f(n)=n.
ｆ０ ＜－ ｆｕｎｃｔｉｏｎ（ｎ）｛
ｙ＝ａｓ．ｃｈａｒａｃｔｅｒ（ｎ）
ｎｃ＝ｎｃｈａｒ（ｙ）
ｚ＝ＮＵＬＬ
ｆｏｒ（ｉ ｉｎ １：ｎｃ） ｚ［ｉ］ ＜－ ｓｕｂｓｔｒ（ｙ，ｉ，ｉ）＝＝１
ｒｅｔｕｒｎ（ｓｕｍ（ｚ））
｝
ｆ ＜－ ｆｕｎｃｔｉｏｎ（ｎ）｛
ｒｅ＝ｎｕｍｅｒｉｃ（）
ｒｅ［１］＝１
ｆｏｒ（ｉ ｉｎ ２：ｎ）｛
ｒｅ［ｉ］＝ｒｅ［ｉ－１］＋ｆ０（ｉ）
｝
ｒｅｔｕｒｎ（ｒｅ［ｎ］）
｝
ｇ ＜－ ｆｕｎｃｔｉｏｎ（ｎ） ｎ － ｆ（ｎ）
ｉ＝１９９９７９
ｗｈｉｌｅ（ｇ（ｉ）！＝０）｛
ｉ＜－ｉ＋１
｝
ｉ
631:卵の名無しさん
18/07/23 22:11:59.00 B1Q5jsGg.net
6!=(5+1)5!=5*5!+5!
=5*5!+(4+1)*4!
=5*5!+4*4!+4!
=5*5!+4*4!+(3+1)*3!
=5*5!+4*4!+3*3!+3!
=5*5!+4*4!+3*3!+(2+1)*2!
=5*5!+4*4!+3*3!+2*2!+2!
=5*5!+4*4!+3*3!+2*2!+1*1!+1
factorial(n+1)= sum(factoral(n:1)*(n:1))+1
632:卵の名無しさん
18/07/24 09:34:20.21 RZUDxWdB.net
Ｂ，Ｃ，Ｄ，Ｅ，Ｆが０～９の数字（同じ数字であってもよい）で
６！＊Ｂ＋５！＊Ｃ＋４！＊Ｄ＋３！＊Ｅ＋２！＊Ｆ＋１！＊Ｇ＝５５５５
が成立するときＢ＋Ｃ＋Ｄ＋Ｅ＋Ｆ＋Ｇの最小となるＢ～Ｇの組み合わせを求めよ。
n=5555
f <- function(B,C,D,E,F,G) sum(factorial(6:1)*c(B,C,D,E,F,G))-n
max=n%/%factorial(6)
gr=expand.grid(0:max,0:6,0:5,0:4,0:3,0:2)
ret=mapply(f, gr[,1],gr[,2],gr[,3],gr[,4],gr[,5],gr[,6])
bg=apply(gr[which(ret==0),],1,sum)
(min.sum=bg[which.min(bg)])
gr[names(min.sum),]
n%/%factorial(6)
n%%factorial(6)%/%factorial(5)
n%%factorial(6)%%factorial(5)%/%factorial(4)
n%%factorial(6)%%factorial(5)%%factorial(4)%/%factorial(3)
n%%factorial(6)%%factorial(5)%%factorial(4)%%factorial(3)%/%factorial(2)
n%%factorial(6)%%factorial(5)%%factorial(4)%%factorial(3)%%factorial(2)%/%factorial(1)
633:卵の名無しさん
18/07/24 16:06:02.90 Cp/BkN3V.net
#include<stdio.h>
#include<string.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;
}
main( int argc, char *argv[] ){
char N = '0'+ atoi(argv[1]);
int a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t;
for(a = 0; a < 10; a++){
for(b = 0; b < 10; b++){
for(c = 0; c < 10; c++){
for(d = 0; d < 10; d++){
for(e = 0; e < 10; e++){
for(f = 0; f < 10; f++){
for(g = 0; g < 10; g++){
for(h = 0; h < 10; h++){
for(i = 0; i < 10; i++){
for(j = 0; j < 10; j++){
for(k = 0; k < 10; k++){
for(l = 0; l < 10; l++){
634:卵の名無しさん
18/07/24 16:07:14.87 Cp/BkN3V.net
for(m = 0; m < 10; m++){
for(n = 0; n < 10; n++){
for(o = 0; o < 10; o++){
for(p = 0; p < 10; p++){
for(q = 0; q < 10; q++){
for(r = 0; r < 10; r++){
for(s = 0; s < 10; s++){
for(t = 0; t < 10; t++){
int num[]={a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t};
qsort(num,20,sizeof(int),compare_int);
char NUM[20];
int i;
for(i=0;i<20;i++){
NUM[i]='0'+ num[i];
}
char N4[4]={N,N,N,N};
char N5[5]={N,N,N,N,N};
if(strstr(NUM,N4)!=NULL & strstr(NUM,N5)==NULL &
(100*a+10*b+c)*f == 100*g+10*h+i &
(100*a+10*b+c)*e == 100*j+10*k+l &
(100*a+10*b+c)*d == 100*m+10*n+o &
(100*m+10*n+o)*100 + (100*j+10*k+l)*10 + 100*g+10*h+i == p*10000+q*1000+r*100+s*10+t
){
printf("abc def ghi jkl mno pqrst = %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d%d%d\n",
a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t);
}
}}}}} }}}}} }}}}} }}}}}
}
635:卵の名無しさん
18/07/28 06:54:00.42 Y3H/zTDn.net
options(scipen = 10)
fibon <- function(N){
f=numeric(N)
f[1]=1
f[2]=1
if(N>2){for(i in 3:N){
f[i]=f[i-2]+f[i-1]
}}
return(f[N])
}
636:Vectorize(fibon)(1:15) fibon(50) fibon(100)# wrong!
637:卵の名無しさん
18/07/31 22:47:13.76 SHi6FXfD.net
c([],L,L).
c([X|L1],L2,[X|L3]) :- c(L1,L2,L3).
638:卵の名無しさん
18/08/01 09:03:07.74 PtTCNKBk.net
?- L= [1,2,3,4,5,6,7,8,9],c([_,_,_],X,L),c(Y,[_,_,_],X).
L = [1,2,3,4,5,6,7,8,9],
X = [4,5,6,7,8,9],
Y = [4,5,6] ;
No
639:卵の名無しさん
18/08/01 21:24:01.69 nhdpq/Mi.net
last([Last],Last).
last([_|Rest],Last) :- last(Rest,Last).
% last([1,2,3,4,5],Last).
640:卵の名無しさん
18/08/04 20:38:03.88 dITWU8BY.net
日本人の血液型はＡ,Ｏ,Ｂ,ＡＢの比率が４：３：２：１であるという。
それぞれの血液型の人を最低でも各々４、３、２、１人集めるためには必要な人数の期待値はいくらか？
641:卵の名無しさん
18/08/04 21:45:28.37 Tkj7u78K.net
>>597
blood.type <- function(p,need){
n=length(p)
ABO=NULL
enough <- function(x){
pool=numeric(n)
for(i in 1:n) pool[i]=sum(ABO==i)
all(pool >= need)
}
while(!enough(ABO)){
ABO=append(sample(1:n,1,p=p),ABO)
}
return(length(ABO))
}
p=4:1
need=c(4,3,2,1)
re=replicate(1e5,blood.type(p,need))
mean(re)
642:卵の名無しさん
18/08/04 21:52:26.83 Tkj7u78K.net
blood.samples <- function(p=c(4,3,2,1),need=c(1,1,1,1),k=1e5){
mean(replicate(k,blood.type(p,need)))
}
blood.samples()
blood.samples(need=4:1,k=1e4)
643:卵の名無しさん
18/08/05 07:20:33.69 aHWo4FJr.net
ｓｕｂｓｅｔ（［］，［］）．
ｓｕｂｓｅｔ（［Ｆｉｒｓｔ｜Ｒｅｓｔ］，［Ｆｉｒｓｔ｜Ｓｕｂ１］） ：－ ｓｕｂｓｅｔ（Ｒｅｓｔ，Ｓｕｂ１）．
ｓｕｂｓｅｔ（［＿｜Ｒｅｓｔ］，Ｓｕｂ２） ：－ ｓｕｂｓｅｔ（Ｒｅｓｔ，Ｓｕｂ２）．
？－ ｓｕｂｓｅｔ（［ド底辺，特殊，シリツ医大，イカサマ入試，裏口，馬鹿］，＿恥）， ｗｒｉｔｅｌｎ（＿恥），ｆａｉｌ．
644:卵の名無しさん
18/08/05 18:08:20.35 cLIwu37J.net
URLﾘﾝｸ(i.imgur.com)
645:卵の名無しさん
18/08/05 20:30:16.90 aHWo4FJr.net
#URLﾘﾝｸ(i.imgur.com)
N=5
library(gtools)
perm=permutations(n=2,r=N,v=c(1,5), repeats.allowed = TRUE)
move=t(apply(perm,1,cumsum))
p0=0
P=(p0+move)%%6
q0=2
Q=(q0+move)%%6
is.fe <- function(x,y){ # is first encounter after N sec?
all(x[-N]!=y[-N]) & x[N]==y[N]
}
re=NULL
for(i in 1:nrow(P)){
for(j in 1:nrow(Q)){
re=append(re,is.fe(P[i,],Q[j,]))
}
}
(Ans=mean(re)) ; 3^(N-1)/(2^N)^2
cat(Ans*(2^N)^2,'/',(2^N)^2)
646:卵の名無しさん
18/08/06 07:05:55.55 29wsmNi+.net
is.pe <- function(x,y,M) { # encounter M or more than M times within N sec.
sum(x==y) >= M
}
M=2
re=NULL
for(i in 1:nrow(P)){
for(j in 1:nrow(Q)){
re=append(re,is.pe(P[i,],Q[j,],M))
}
}
(Ans=mean(re))
cat(Ans*(2^N)^2,'/',(2^N)^2)
647:卵の名無しさん
18/08/08 19:27:51.20 OkpyQ+1n.net
URLﾘﾝｸ(ailaby.com)
648:卵の名無しさん
18/08/08 19:31:00.52 OkpyQ+1n.net
hanoi <- function(n,from='A',via='B',to='C'){
if(n >= 1){
Recall(n-1,from,to,via)
cat('move',n,'from',from,' to ',to,'\n')
Recall(n-1,via,from,to)
}
}
hanoi(4)
649:卵の名無しさん
18/08/09 05:35:20.13 mD+P4tQf.net
掛け算を再帰関数で定義する。
VAT <- function(n,m=1.08){
if(n==0) 0
else Recall(n-1) + m
}
VAT(250)
650:卵の名無しさん
18/08/09 06:51:20.87 mD+P4tQf.net
総和を再帰関数で書く
sup <- function(v){
if (length(v)==0) 0
else v[1] + Recall(v[-1])
}
sup(1:10)
651:卵の名無しさん
18/08/09 10:57:33.87 JutZ+/A9.net
# Ackermann
> A <- function(m,n){
+ if(m==0) return(n+1)
+ if(n==0) return(Recall(m-1,1))
+ else return(Recall(m-1,Recall(m,n-1)))
+ }
>
> A(2,4)
[1] 11
> A(3,4)
[1] 125
> A(4,1)
Error: evaluation nested too deeply: infinite recursion / options(expressions=)?
>
652:卵の名無しさん
18/08/09 18:49:05.66 JutZ+/A9.net
累積和
# cumsum with for-loop
cumsumL <- function(v){
n=length(v)
re=numeric(n)
re[1]=v[1]
for(i in 1:(n-1)) re[i+1]=re[i]+v[i+1]
re
}
# cumsum with recursive call
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)
}
}
}
653:卵の名無しさん
18/08/09 20:12:53.52 JutZ+/A9.net
# 10進法をN進法でdigit桁表示する
dec2n <- function(num, N = 2, digit = 0){ # decimal to 0,1,..,n-1 vector
if(num==0 & digit==0) return(0)
if(num <= 0 & digit <= 0) return()
else{
return(append(Recall(num%/%N, N ,digit-1), num%%N))
}
}
> dec2n(0)
[1] 0
> dec2n(11,digit=5)
[1] 0 0 1 0 1 1
> dec2n(9,N=5,digit=3)
[1] 0 0 1 4
> dec2n(1000,N=16)
[1] 3 14 8
654:卵の名無しさん
18/08/09 21:28:39.07 JutZ+/A9.net
>>610 (degugged)
dec2n <- function(num, N = 2, digit = 1){ # decimal to 0,1,..,n-1 vector
if(num <= 0 & digit <= 0) return()
else{
return(append(Recall(num%/%N, N ,digit-1), num%%N))
}
}
dec2n(0)
dec2n(11,digit=5)
dec2n(9,N=5,digit=3)
dec2n(1000,N=16)
#
dec2hex <- function(x){ # decimal to hexa
hex=c(0:9,letters[1:6])
n=length(x)
re=numeric(n)
for(i in 1:n){
if(x[i]==0) re[i]='0'
else re[i]=hex[x[i]+1]
}
cat(re,'\n')
}
dec2hex(c(1,0,15))
dec2hex(dec2n(1000,16))
dec2sexa <- function(x) dec2hex(dec2n(x,16))
655:卵の名無しさん
18/08/14 14:27:35.01 Z2jjlChF.net
draft
X <- function(n,red=10,white=90){
rw=red+white
total=martix(0,nrow=n,ncol=2^n)
p=martix(0,nrow=n,ncol=2^n)
total[1,1:2]=c(rw-1,rw+1)
p[1,1:2]=c(white/rw,red/rw)
if(n > 1){
for(i in 1:n){
li=2^i
total[i+1,1:(2*li)]=c(total[i,1:li]-1, total[i,1:li]+1)
p[i+1,1:2*li)]=c(p[i,1:li]*(total[i,1:li]-red)/total[i,1:li], p[i,1:li]*(red)/total[i,1:li])
}
}
return(sum(p[n,]*total[n,]))
}
656:卵の名無しさん
18/08/15 08:08:12.45 SmB+loM1.net
rm(list=ls())
X = function(n,r=10,w=90){ #n:試行数 r: 赤玉数 w:白玉数
rw=r+w # 試行前総玉数
J=rw+n # 総玉数の上限
# s[i,j] i回試行後に総数がj個である確率の行列
s=matrix(0,nrow=n,ncol=J)
s[1,rw-1]=w/rw ; s[1,rw+1]=r/rw # 1回試行後
if(n > 1){
for(i in 2:n){
for(j in r:J){ # jはr未満にはならない
# if(j==1) s[i,j] = s[i-1,j+1]*(j+1-r)/(j+1)
if(j==J) s[i,j] = s[i-1,j-1] * r/(j-1)
else s[i,j] = s[i-1,j-1] * r/(j-1) + s[i-1,j+1]*(j+1-r)/(j+1)
}
}
}
total=sum((r:J)*s[n,r:J])# n回試行後総数の期待値
white=total-r
return(c(total=total,white=white))
}
> vX=Vectorize(X)
> vX(1000:1010)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
total 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5
white 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5
657:卵の名無しさん
18/08/15 12:17:01.96 eDY+Gzxn.net
# シオマネキ
# URLﾘﾝｸ(i.imgur.com)
H=(sqrt(5+2*sqrt(5)))/2
A=0.5+0i
B=A+r
shiomaneki <- function(x,y){
xy=x+y*1i
(abs(xy-A)+abs(xy-B)+x)*2 + 1
}
x=seq(0,Re(B),length=100)
y=seq(0,H,length=100)
z=outer(x,y,shiomaneki)
contour(x,y,z,nlevels=20,lwd=2,bty='l')
x11() ;persp(x,y,z, theta=35, lty=3,col="lightblue",xlab='x',ylab='y',zlab='z',
ticktype='detailed',shade=0.75,phi=30,ltheta=-10,border=TRUE)
rgl::persp3d(x,y,z,col='blue')
sio <- function(xy) shiomaneki(xy[1],xy[2])
optim(c(0.4,0.4),sio,method='Nelder-Mead')
658:卵の名無しさん
18/08/15 12:17:36.49 eDY+Gzxn.net
# カブトガニ
# URLﾘﾝｸ(i.imgur.com)
H=(sqrt(5+2*sqrt(5)))/2
A=0.5+0i
B=A+r
C=0+H*1i
abs(C-B)
kabutogani <- function(y1,y2){
K=0+y1*1i
G=0.5+y2*1i
abs(K-C) + 2*abs(K-G) + 2*abs(G-A) +2*abs(G-B)
}
y1=seq(0,H,length=100)
y2=seq(0,H,length=100)
z=outer(y1,y2,kabutogani)
contour(y1,y2,
659:z,nlevel=100,lwd=2,xlim=c(0.5,1.5),ylim=c(0,1.2),bty='l') kab <- function(y1y2) kabutogani(y1y2[1],y1y2[2]) optim(c(1,0.5),kab, method = 'Nelder-Mead')
660:卵の名無しさん
18/08/15 15:14:03.13 eDY+Gzxn.net
DOP <- function(n,print=FALSE){ # diagonal length of regular plygon of size length 1
q=2*pi/n
r=cos(q)+1i*sin(q)
p=numeric(n+1)
for(i in 1:(n+1)) p[i]=r^(i-1)
D=NULL
if(n>3){
for(j in 3:(2+ceiling(n-3)/2)){
D= append(D,abs(p[1]-p[j]))
}
}
if(print){
plot(p,type='l',bty='l',axes=FALSE,ann=FALSE,lwd=2)
for(i in 3:(n-1)) segments(1,0,Re(p[i]),Im(p[i]),col='gray')
}
return(D)
}
DOP(17,p=T)
661:卵の名無しさん
18/08/15 22:04:58.96 SmB+loM1.net
DOP <- function(n,print=FALSE){ # diagonal length of regular plygon of size length 1
q=2*pi/n
r=cos(q)+1i*sin(q)
p=numeric(n+1)
for(i in 1:(n+1)) p[i]=r^(i-1)
D=NULL
if(n>3){
for(j in 3:(2+ceiling(n-3)/2)){
D= append(D,abs(p[1]-p[j])/abs(p[1]-p[2]))
}
}
if(print){
plot(p,type='l',bty='l',axes=FALSE,ann=FALSE,lwd=2)
for(i in 3:(n-1)) segments(1,0,Re(p[i]),Im(p[i]),col='gray')
}
return(D)
}
DOP(17,p=T)
662:卵の名無しさん
18/08/15 23:51:51.61 eDY+Gzxn.net
> (sin(2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - (cos(4*2*pi/7)*sin(2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - sin(4*2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) + (cos(2*pi/7)* sin(4*2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2))
[1] 2.190643
663:卵の名無しさん
18/08/16 22:05:33.14 8kKe3yXf.net
rm(list=ls())
graphics.off()
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,lwd=1,asp=1)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),1:n,cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
664:卵の名無しさん
18/08/17 09:11:30.66 qXmvzl8y.net
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),paste('p',1:n),cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),col=2,lwd=2,...)
pt <- function(x,y,...) text(Re(x),Im(x), y,...)
kabutogani3 <- function(xl,yl,xc,yc,xr,yr){
L=xl+yl*1i
C=xc+yc*1i
R=xr+yr*1i
abs(p[3]-C)+abs(C-R)+abs(C-L)+abs(R-p[2])+abs(L-p[4])+abs(R-p[1])+abs(L-p[5])
}
kabu3 <- function(par){
kabutogani3(par[1],par[2],par[3],par[4],par[5],par[6])
}
p=ngon(5,axis=T,col='lightblue',lwd=2)
opt=optim(runif(6),kabu3,method='Nelder-Mead')
kabu3(opt$par)
(par=opt$par)
L=par[1]+par[2]*1i
C=par[3]+par[4]*1i
R=par[5]+par[6]*1i
pt(C,'C') ; pt(L,'L') ; pt(R,'R')
seg(p[1],R);seg(p[2],R);seg(p[3],C);seg(p[4],L);seg(p[5],L);seg(L,C);seg(R,C)
665:卵の名無しさん
18/08/17 09:24:47.60 qXmvzl8y.net
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),paste('p',1:n),cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),col=2,lwd=2,...)
pt <- function(x,y,...) text(Re(x),Im(x), y,...)
kabutogani3 <- function(xl,yl,xc,yc,xr,yr){
L=xl+yl*1i
C=xc+yc*1i
R=xr+yr*1i
abs(p[3]-C)+abs(C-R)+abs(C-L)+abs(R-p[2])+abs(L-p[4])+abs(R-p[1])+abs(L-p[5])
}
kabu3 <- function(par){
kabutogani3(par[1],par[2],par[3],par[4],par[5],par[6])
}
p=ngon(5,axis=T,col='lightblue',lwd=2)
opt=optim(runif(6),kabu3,method='CG')
kabu3(opt$par)
(par=opt$par)
L=par[1]+par[2]*1i
C=par[3]+par[4]*1i
R=par[5]+par[6]*1i
pt(C,'C') ; pt(L,'L') ; pt(R,'R')
seg(p[1],R);seg(p[2],R);seg(p[3],C);seg(p[4],L);seg(p[5],L);seg(L,C);seg(R,C)
666:卵の名無しさん
18/08/17 17:56:05.86 qXmvzl8y.net
# how many ways of allocating 5 rooms to 6 people without vacancy?
# allocated to 1 room (4 vacant)
a1=choose(5,1)*1^6 ; a1
# allocated to 2 rooms (3 vacant)
a2=choose(5,2)*(2^6-2) ; a2
# allocated to 3 rooms (2 vacant)
a3=choose(5,3)*( 3^6-choose(3,2)*(2^6-2)-3 ) ; a3
# allocated to 4 rooms (1 vacant)
a4=choose(5,4)*( 4^6 - choose(4,3)*(3^6-choose(3,2)*(2^6-2)-3) - choose(4,2)*(2^6-2)-4 ) ; a4
5^6 - a1 - a2 - a3 - a4
667:卵の名無しさん
18/08/17 18:38:19.94 qXmvzl8y.net
dec2n <- function(num, N = 2, digit = 1){ # decimal to 0,1,..,n-1 vector
if(num <= 0 & digit <= 0) return()
else{
return(append(dec2n(num%/%N, N ,digit-1), num%%N))
}
}
room.allocation <- function(n,r){ # allocate n people to r rooms without vacancy
max=r^n
counter=0
x=0
while(x < max){
if(length(unique(dec2n(x,r,n))) == r) counter = counter+1
x= x + 1
}
return(counter)
}
668:卵の名無しさん
18/08/18 16:03:38.91 Z1UCnKoz.net
kagamimochi = function(b, h){
r = (b^2 + h^2)/(2*h)
V = (2/3*r^3 - r^2*(r-h) + (r-h)^3/3)*pi
return(c(radius=r,Volume=V))
}
669:卵の名無しさん
18/08/18 17:38:21.87 Z1UCnKoz.net
kagamimochi = function(b, h){
r = (b^2 + h^2)/(2*h)
if(b > h) V = (2/3*r^3 - r^2*(r-h) + (r-h)^3/3)*pi
else V = (r^2*(h-r) - 1/3*(h-r)^3 + 2*r^3/3)*pi
return(c(radius=r,Volume=V))
}
670:卵の名無しさん
18/08/18 20:09:40.32 TnH4H/8z.net
draft
# how many ways of allocating 6 rooms to n people without vacancy?
# allocated to 1 room
a1=choose(6,1)*1^n
# allocated to 2 rooms
a2=choose(6,2)*(2^n-2)
# allocated to 3 rooms
a3=choose(6,3)*( 3^n-choose(3,2)*(2^n-2)-3 )
# allocated to 4 rooms
a4=choose(6,4)*( 4^n - choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4 )
# allocated to 5 rooms
a5=choose(6,5)*( 5^n - choose(5,4)*(4^n-choose(4,3)*(3^n-
choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4)-choose (5,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose (5,2)*(2^n -2) -5
6^n - a1 - a2 - a3 - a4 - a5
671:卵の名無しさん
18/08/18 20:22:13.81 Z1UCnKoz.net
# how many ways of allocating 6 rooms to n people without vacancy?
library(Rmpfr)
six_rooms <- function(x){
if(x[1]<10) n=x else n=mpfr(x,100)
# allocated to 1 room
a1=choose(6,1)*1^n
# allocated to 2 rooms
a2=choose(6,2)*(2^n-2)
# allocated to 3 rooms
a3=choose(6,3)*( 3^n-choose(3,2)*(2^n-2)-3 )
# allocated to 4 rooms
a4=choose(6,4)*( 4^n - choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4 )
# allocated to 5 rooms
a5=choose(6,5)*(5^n-choose(5,4)*(4^n-choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3)
-choose(4,2)*(2^n-2)-4)-choose(5,3)*(3^n-choose(3,2)*(2^n-2)-3)-choose(5,2)*(2^n -2) -5)
6^n - a1 - a2 - a3 - a4 - a5
}
672:卵の名無しさん
18/08/19 00:27:23.03 gaGSkZ47.net
allocate.rooms <- function(m,n){ # m:rooms n:people
if(m==n) return(factorial(m))
else if(m==1) return(1)
else m*Recall(m,n-1) + m*Recall(m-1,n-1)
}
#include<stdio.h>
long factorial(long n) {
long re = 1;
long k;
for(k=1;k <=n;k++) {re *= k;}
return re;
}
long rooms(int m, int n){
if(m==n) { return factorial(m);}
else if(m==1){ return 1;}
else{
return m * rooms(m,n-1) + m * rooms(m-1,n-1);
}
}
void main( int argc, char *argv[] ){
int m,n;
long ways;
m=atoi(argv[1]);
n=atoi(argv[2]);
ways=rooms(m,n);
printf("%d\n",ways);
}
673:卵の名無しさん
18/08/19 15:49:26.54 gaGSkZ47.net
楕円体 x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
高さ h のロケットおっぱいの体積
RocketPi <= function(a,b,c,h) 2/3*pi*a*b*c - 1/3*pi*b*c*(3*a^2-(a-h)^2)*(a-h)/a
674:卵の名無しさん
18/08/19 15:56:15.35 gaGSkZ47.net
楕円体 x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
RocketPi <= function(a,β,γ,h) # a:楕円体の長軸長, β,γ:ロケットおっぱいの楕円底面の軸長,ｈ：ロケットおっぱいの高さ
b=(a-h)*β/sqrt(a^2-(a-h)^2)
c=(a-h)*γ/sqrt(a^2-(a-h)^2)
2/3*pi*a*b*c - 1/3*pi*b*c*(3*a^2-(a-h)^2)*(a-h)/a
}
675:卵の名無しさん
18/08/19 17:09:55.54 gaGSkZ47.net
# x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
(1+sqrt(5))/2
f <- function(a=10,b=10,c=10){
xy2z <- function(x,y) c*sqrt(a^2*b^2-a^2*y^2-b^2*x^2)/(a*b)
x=seq(-a,a,le=50)
y=seq(-r*a,r*a,le=50)
z=outer(x,y,xy2z)
contour(x,y,z)
persp(x,y,z,
theta=35, lty=3,col="pink",xlab='x',ylab='y',zlab='z',ylim=c(-2*a,2*a),
ticktype='detailed',shade=0.75,phi=30,ltheta=-10,border=TRUE)
rgl::persp3d(x,y,z,col='pink')
}
f()
f(10,20,25)
676:卵の名無しさん
18/08/19 23:40:10.33 gaGSkZ47.net
seg <- function(a,b,...){
segments(Re(a),Im(a),Re(b),Im(b),...)}
pt <- function(x,y=NULL,...){
text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)}
# solve α^2/a^2 + (b+β)^2/b^2=1 for a
α=3
β=-1
b=1
# α^2*b^2/(-β*(2*b+β))
(a = α*b/sqrt(-β*(2*b+β)))
(q=b+β)
f= function(x,y) x^2/a^2 + (y-q)^2/b^2 # = 1
x=seq(-5,5,length=100)
y=seq(-5,5,length=100)
z=outer(x,y,f)
contour(x,y,z,level=1,bty='l',xlim=c(-10,10),ylim=c(-2,10),
drawlabels=FALSE,col=sample(colours(),1),axes=FALSE,lwd=2)
pt(α,'●') ; pt(-α,'●') ; pt(β*1i,'●')
seg(-10,10,lty=3) ; seg(-4i,10i,lty=3)
pt(0+q*1i,'.')
eclipse <- function(b) {
# α^2*b^2/(-β*(2*b+β))
(a = α*b/sqrt(-β*(2*b+β)))
(q=b+β)
f= function(x,y) x^2/a^2 + (y-q)^2/b^2 # = 1
x=seq(-10,10,length=100)
y=seq(-10,10,length=100)
z=outer(x,y,f)
contour(x,y,z,level=1,drawlabels=FALSE,add=TRUE,col=sample(colours(),1),lwd=2)
pt(0+q*1i,'.')
}
for(i in 2:10) eclipse(i)