17/10/09 08:02:11.41 Wked2yci.net
# 10分の1で10万円が当たるクジAと1万分の1で8000万円が当たるクジの方Bある
# クジ1本の値段は前者は2万円,後者は4万円とする.
# 購入金額100万円としてどのように買うとき賞金の期待値が最大になるか？
pa=0.1 #Aの当選率
Wa=10 #Aの賞金
ca=2 #Aのコスト
pb=1/10000
Wb=8000
cb=4
Cash=100 #購入額
# a,bは各々A,Bの購入数 2*a+4*b=100 a=(Cash-4*b)/2
award <- function(b){
a=(100-4*b)/2
sum(rbinom(a, 1, pa) * Wa) + sum(rbinom(b, 1, pb) * Wb)
}
bb=0:floor(Cash/cb)
AWD <- function(){
Award=sapply(bb,award)
c(which.max(Award)-1,max(Award))
}
k=10^6 ; res=replicate(k,AWD())
Max.B=res[1,]
hist(Max.B,freq=FALSE,breaks=20,col='lightblue')
lines(density(Max.B),lwd=1)
summary(Max.B)
Mode(Max.B)
Max.AWD=res[2,]
hist(Max.AWD,col='wheat',breaks=100)
summary(Max.AWD)
Mode(Max.AWD)
136:卵の名無しさん
17/10/09 11:07:15.42 Wked2yci.net<
137:> Rstudioのdplyrとtidyr の Cheetsheet https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf のデータは EDAWR というパッケージにあるらしいのだが、R3.4.2だと Warning in install.packages : package ‘EDAWR’ is not available (for R version 3.4.2) と警告がでた。 install.packages('devtools') devtools::install_github("rstudio/EDAWR") で警告を回避してインストールできた。 > cases country 2011 2012 2013 1 FR 7000 6900 7000 2 DE 5800 6000 6200 3 US 15000 14000 13000 > tidyr::gather(cases,'year','n',2:4) country year n 1 FR 2011 7000 2 DE 2011 5800 3 US 2011 15000 4 FR 2012 6900 5 DE 2012 6000 6 US 2012 14000 7 FR 2013 7000 8 DE 2013 6200 9 US 2013 13000 と動作してくれた。
138:卵の名無しさん
17/10/09 14:34:35.08 Wked2yci.net
library(tidyr)
library(EDAWR)
cases
cases %>% gather('year','n',2:4)
cases %>% gather('year','n',2:4) %>% spread(year,n) %>% arrange(country)
cases %>% arrange(country)
pollution
pollution %>% spread(size, amount)
pollution %>% spread(size, amount) %>% gather('size','amount',2:3) %>% arrange(desc(city))
pollution
library(reshape2)
cases
cases %>% melt(value.name ='cases')
pollution
pollution %>% acast(city ~ size)
139:innuendo
17/10/09 15:13:49.58 N8CcYMr4.net
>>127
URLﾘﾝｸ(imagizer.imageshack.com)
140:卵の名無しさん
17/10/09 15:39:38.20 Wked2yci.net
>>135
URLﾘﾝｸ(i.imgur.com)
141:卵の名無しさん
17/10/09 15:40:36.46 Wked2yci.net
> summary(g13)
Min. 1st Qu. Median Mean 3rd Qu. Max.
8992 9866 9933 9903 9972 10000
> which.max(table(g13))
10000
714
142:卵の名無しさん
17/10/09 19:26:11.72 Wked2yci.net
pa=0.1 ; pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
k=10^3 # シミュレーション回数
# Aをi買ったときの賞金総額
a <- function(i){sum(rbinom(i, 1, pa) * Wa) + sum(rbinom(100-i, 1, pb) * Wb)}
which.max2 <- function(x){which(x==max(x))-1}
# 最大賞金総額になるAの購入数（複数可）の配列を返す
b <- function(){
b1=NULL
b1=c(b1,which.max2
143:(sapply(0:n, a))) return(b1) } b2=NULL for(j in 1:k){ b2=c(b2,b()) } hist(b2,freq=FALSE,breaks=20,col=sample(colors(),2),ann=FALSE) lines(density(b2),lwd=1) summary(b2) which.max(table(b2))
144:卵の名無しさん
17/10/10 06:51:19.25 wISdLCFn.net
pa=0.1 ; pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
aa=0:n # 宝くじAを買う候補
Get <- function(a){sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb} #賞金総額
foo <- function(k=10^3){
which.max2 <- function(x){which(x==max(x))-1}
MAX=NULL
MAX.AWD=NULL
for(i in 1:k){
tmp=sapply(aa,Get)
indx=which.max2(tmp)
MAX=c(MAX,indx) # 最大賞金を獲得したときに買ったAの本数の配列
MAX.AWD=c(MAX.AWD,tmp[indx+1]) #最大賞金が複数回あるときは重複して数える
}
indx=which.max2(MAX.AWD)+1
plot(MAX,MAX.AWD,xlim=c(0,n),ylim=c(n*pa*Wa,8*n*pa*Wa),
xlab='A_lots',ylab='Award',col=rgb(.1,.1,.1,.5))
points(MAX[indx],MAX.AWD[indx],pch=19)
luckyA=MAX[indx]
lucky.AWD=MAX.AWD[indx]
cbind(luckyA,lucky.AWD)
}
foo()
145:卵の名無しさん
17/10/10 16:36:50.99 g2TXzqHS.net
> '+'(1,2)
[1] 3
> '*'(2,3)
[1] 6
> a=1:10
> '['(a,7)
[1] 7
'[' って関数なんだな。
146:卵の名無しさん
17/10/10 21:23:23.64 VF4JdnZT.net
>>129
宝くじAは1000本中300本が当たりで当たりは 100万円の賞金、
宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、
とAの期待値＝Bの期待値の3倍のとき
URLﾘﾝｸ(i.imgur.com)
147:卵の名無しさん
17/10/11 06:48:41.73 8dgUXTG5.net
>>139
続きのコードが書き込めないので画像でアップ
URLﾘﾝｸ(imagizer.imageshack.com)
148:卵の名無しさん
17/10/11 07:14:30.46 8dgUXTG5.net
# pa=0.1
pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
foo <- function(pa=0.1,k=10^3){
# k=10^3 # シミュレーション回数
# Aをi買ったときの賞金総額
a <- function(i){sum(rbinom(i, 1, pa) * Wa) + sum(rbinom(100-i, 1, pb) * Wb)}
which.max2 <- function(x){which(x==max(x))-1}
# 最大賞金総額になるAの購入数（複数可）の配列を返す
b <- function(){
b1=NULL
b1=c(b1,which.max2(sapply(0:n, a)))
return(b1)
}
b2=NULL
for(j in 1:k){
b2=c(b2,b())
}
hist(b2,freq=FALSE,breaks=20,col=sample(colors(),1),xlab='',ylab='',yaxt='n',
main=paste('pa =',pa))
lines(density(b2),lwd=1)
print(summary(b2))
print(c(Mode=as.numeric(names(which.max(table(b2)))))) # モード値
print(c(length=length(b2))) # 最大値となる数は複数最大値があるのでｎを超える
invisible(b2)
}
par(mfrow=c(2,2))
for(i in c(0.15,0.2,0.25,0.3)) foo(i,500)
149:卵の名無しさん
17/10/11 20:53:31.06 8dgUXTG5.net
確かにマジックだな
##-- now some "magic" to do the 4 regressions in a loop:
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
## or ff[[2]] <- as.name(paste0("y", i))
## ff[[3]] <- as.name(paste0("x", i))
mods[[i]] <- lmi <- lm(ff, data = anscombe)
print(anova(lmi))
}
## See how close they are (numerically!)
sapply(mods, coef)
lapply(mods, function(fm) coef(summary(fm)))
## Now, do what you should have done in the first place: PLOTS
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma = c(0, 0, 2, 0))
for(i in 1:4) {
ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
xlim = c(3, 19), ylim = c(3, 13))
abline(mods[[i]], col = "blue")
}
mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5)
par(op)
150:卵の名無しさん
17/10/11 21:31:48.49 8dgUXTG5.net
[ って関数だったんだ
(d=outer(11:19,11:19,'*'))
d[,3]
apply(d,1,'[',3)
151:卵の名無しさん
17/10/11 23:59:07.55 8dgUXTG5.net
(d=array(1:48, dim=c(6, 4, 2)))
# を，3行2列おきに平均して
# ,,1
# [,1] [,2]
# [1,] 5 17
# [2,] 8 20
# ,,2
# [,1] [,2]
# [1,] 29 41
# [2,] 32 44
#
a=3
b=2
c=1
l=6
m=4
n=2
re=NULL
for(k in 1:(n/c)){
for(j in 1:(m/b)){
for(i in 1:(l/a)){
re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),1:c+c*(k-1)]))
}
}
}
array(re,dim=c(m/b,l/a,n))
152:卵の名無しさん
17/10/12 18:41:56.42 IVfk5sKl.net
fs<-function(){ #スパゲッティ・ソース
d=array(1:(l*m*n),dim=c(l,m,n))
re=NULL
for(k in 1:(n/c)){
for(j in 1:(m/b)){
for(i in 1:(l/a)){
re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),1:c+c*(k-1)]))
}}}
array(re,dim=c(m/b,l/a,n))
}
fe <- function(){ # Expertソース
A <- array(1:(l*m*n), dim=c(l,m,n))
B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c))
f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])}
re=mapply(f, B[,1], B[,2],B[,3])
array(re,dim=c(m%/%b,l%/%a,n%/%c))
}
a=4
b=2
c=1
l=400
m=200
n=2
> system.time(fs())
user system elapsed
2.17 0.00 2.18
> system.time(fe())
user system elapsed
1.14 0.00 1.14
153:卵の名無しさん
17/10/14 04:34:05.51 IFvRV3Gd.net
x=11:19
y=11:19
nx=length(x)
ny=length(y)
re=matrix(rep(NA,nx*ny),ncol=nx)
for(i in 1:nx){
for(j in 1:ny){
re[j,i]=x[i]*y[j]
}
}
re
a=expand.grid(x,y)
b=mapply('*',a[,1],a[,2])
matrix(b,nx)
outer (x,y,'*')
154:卵の名無しさん
17/10/14 06:39:17.90 IFvRV3Gd.net
x=1:999
y=1:999
nx=length(x)
ny=length(y)
fs = function (){
re=matrix(rep(NA,nx*ny),ncol=nx)
for(i in 1:nx){
for(j in 1:ny){
re[j,i]=x[i]*y[j]
}
}
re
}
fe = function (){
a=expand.grid(x,y)
b=mapply('*',a[,1],a[,2])
matrix(b,nx)
}
fo = function (){
outer (x,y,'*')
}
system.time(fs())
system.time(fe())
system.time(fo())
155:卵の名無しさん
17/10/14 19:44:47.68 +G1xzFga.net
##
m=7 ; n=17
Z=matrix(rep(NA,m*n),m,n)
colnames(Z)=paste0('n',0:(n-1))
rownames(Z)=paste0('m',0:(m-1))
f <- function(x) c(x%%m+1,x%%n+1)
for(i in 0:(n*m-1)){
idx=f(i) ; idx
Z[idx[1],idx[2]]=i
}
length(unique(Z))
Z
NZ <- function(x,y,m=7,n=17){
f <- function(x) c(x%%m+1,x%%n+1)
for(i in 0:(n*m-1)){
idx=f(i)
Z[idx[1],idx[2]]=i
}
Z[(x%%m+1),(y%%n+1)]
}
NZ(5,15)
NZ(6,0)
a=expand.grid(0:6,0:16)
res=mapply(NZ,a[,1],a[,2])
Z=matrix(res,m,n)
colnames(Z)=paste0('n',0:(n-1))
rownames(Z)=paste0('m',0:(m-1))
Z
length(unique(Z))
156:卵の名無しさん
17/10/14 19:45:01.23 +G1xzFga.net
where2 <- function(x,m=7,n=17){
modm=paste0('(mod ',m,')')
modn=paste0('(mod ',n,')')
res=c(x%%m,x%%n)
names(res)=c(modm,modn)
return(res)
}
where2(100)
157:卵の名無しさん
17/10/14 19:45:34.01 +G1xzFga.net
## 和で直積
fij <- function(i,j){
ai=where2(i)
aj=where2(j)
aiaj=ai+aj
NZ(aiaj[1],aiaj[2]) == i+j
}
a=expand.grid(0:(m-1),0:(n-1))
mapply(fij,a[,1],a[,2])
## 積で直積
fij <- function(i,j){
ai=where2(i)
aj=where2(j)
aiaj=ai * aj
NZ(aiaj[1],aiaj[2]) == i * j
}
a=expand.grid(0:(m-1),0:(n-1))
mapply(fij,a[,1],a[,2])
158:卵の名無しさん
17/10/14 19:46:23.30 +G1xzFga.net
> mapply(fij,a[,1],a[,2])
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[18] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[35] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[52] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[69] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[86] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[103] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
同型写像確認
159:卵の名無しさん
17/10/14 22:27:16.20 +G1xzFga.net
NZ <- function(x,y,z, l=13,m=15,n=17){
lmn=l*m*n
N=0:(lmn-1)
.M=rbind(nl=N%%l,nm=N%%m,nn=N%%n)
which(apply((.M - c(x,y,z)),2,function(x) sum(x^2))==0)-1
}
l=13 ; m=15 ; n=17
ll=0:(l-1); mm=0:(m-1) ; nn=0:(n-1)
a=expand.grid(ll,mm,nn)
res=mapply(NZ,a[,1],a[,2],a[,3])
length(unique(res)) == l*m*n
summary(res)
plot(sort(res))
hist(res,breaks=l*m*n)
l*m*nNZ <- function(x,y,z, l=13,m=15,n=17){
lmn=l*m*n
N=0:(lmn-1)
.M=rbind(nl=N%%l,nm=N%%m,nn=N%%n)
which(apply((.M - c(x,y,z)),2,function(x) sum(x^2))==0)-1
}
l=13 ; m=15 ; n=17
ll=0:(l-1); mm=0:(m-1) ; nn=0:(n-1)
a=expand.grid(ll,mm,nn)
res=mapply(NZ,a[,1],a[,2],a[,3])
length(unique(res)) == l*m*n
summary(res)
plot(sort(res))
hist(res,breaks=l*m*n)
l*m*n 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)
160:卵の名無しさん
17/10/15 07:15:02.95 hnCh54BK.net
URLﾘﾝｸ(en.wikipedia.org)
を参考にオイラー関数をグラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
phi <- function(n){
nn=0:(n-1)
f=function(x,y) (x*y)%%n
names(nn)=paste0('n',nn)
z=outer(nn,nn,f)
i=which(gmp::gcd(1:(n-1),n)==1)
length(i)
}
#Fourier transform
phi.F <- function(N){
nn=1:N
sum( gmp::gcd(nn,N)*cos(2*pi*nn/N))
}
N=1000
nn=1:N
y=sapply(nn,phi)
plot(nn,y,col=sample(colours()),pch=19,ylab='',xlab='n',main='Euler\'s totient function')
y1=sapply(nn,phi.F)
points(nn,y1,pch='.',cex=1.5)
161:卵の名無しさん
17/10/16 14:46:14.17 SUEY/256.net
(p+1)^(p^(n-1)) ≡1 (mod p^n)
(p+1)^(p^(n-1)) ≡ p^n + 1 (mod p^(n+1))
162:卵の名無しさん
17/10/16 18:24:23.14 SUEY/256.net
# (p+1)^(p^(n-1)) ≡1 (mod p^n)
beki <- function(x,y,p){ # (x^y)%%p
if(y==0) return(1)
if(y==1) return(x%%p)
re=numeric(y)
re[1]=x
for(i in 1:(y-1)) re[i+1] = (x*re[i])%%p
return(re[y])
}
f <- function(p,n){ # (p+1)^(p^(n-1)) ≡1 (mod p^n)
beki(p+1,p^(n-1),p^n)
}
p=(2:9)[gmp::isprime(2:9)!=0]
n=2:9
a=expand.grid(p,n)
> mapply(f,a[,1],a[,2])
> mapply(f,a[,1],a[,2])
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)
163:卵の名無しさん
17/10/16 19:48:05.23 SUEY/256.net
(p+1)^(p^(n-1)) ≡ p^n + 1 (mod p^(n+1))
p:奇素数
164:卵の名無しさん
17/10/18 20:51:15.46 r7jfJx+N.net
# (x^y)%%p
beki <- function(x,y,p){
if(!y) return(x^y)
re=numeric()
re[1]=x%%p
for(i in 1:y) re[i+1] = (x*re[i])%%p
return(re[y])
}
# ｎのmod p での位数を求める
isu <- function(n,p){
re=numeric()
for(i in 1:(p-1)){
re[i]=beki(n,i,p)
}
which.min(re)
}
165:卵の名無しさん
17/10/19 15:04:12.35 lylsDJ3i.net
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
e =c(1,2,3)
d =c(3,1,2)
d2=c(2,3,1)
t =c(1,3,2)
td=c(3,2,1)
td2=c(2,1,3)
(.M=rbind(e,d,d2,t,td,td2))
(names=rownames(.M))
equal <- function(x,y) sum((x-y)^2)==0
f <- function(x){ #
re=NULL
for(i in 1:nrow(.M)){
if(equal(.M[i,],x)) return(i) # i for names[i]
}
}
g <- function(x,y){
f(tikan2(x,y)) # first y, then x makes names[i]
}
166:卵の名無しさん
17/10/19 15:04:48.68 lylsDJ3i.net
.L = list(e,d,d2,t,td,td2)
names(.L)=names
.L
a=expand.grid(.L,.L)
(.a=mapply(g,a[,1],a[,2]))
(.M0=matrix(mapply(g,a[,1],a[,2]),nrow(.M)))
.M1=matrix(names[.a], nrow(.M))
rownames(.M1)=names
colnames(.M1)=names
> print(.M1, quote=FALSE)
e d d2 t td td2
e e d d2 t td td2
d d d2 e td2 t td
d2 d2 e d td td2 t
t t td td2 e d d2
td td td2 t d2 e d
td2 td2 t td d d2 e
167:卵の名無しさん
17/10/19 16:00:22.26 lylsDJ3i.net
結合側の検証
k <- function(x,y,z) equal(tikan2(x,tikan2(y,z)), tikan2(tikan2(x,y),z))
(a=expand.grid(.L,.L,.L))
b=mapply(k,a[,1],a[,2],a[,3])
> summary(b)
Mode TRUE
logical 216
168:卵の名無しさん
17/10/19 19:20:37.10 pl4P+6GE.net
>>162
結合則
169:卵の名無しさん
17/10/19 19:26:48.29 ZZphLwM2.net
演算が閉じている、単位元が存在する、逆元が存在する、の確認は容易だけど
結合則の確認は手計算だと手間がかかる。
プログラムを組んで確認できて納得。
170:卵の名無しさん
17/10/20 13:57:10.48 hxLmqOB4.net
紛らわしいなぁ
群論における「位数」とは、二つの異なる定義があり、それぞれ、
①有限群Gの位数＝有限群Gの元の個数、
②有限群Gの元aの位数＝有限群Gの元aにおいて、a^m=e となる最小の正の整数、
となります。
171:卵の名無しさん
17/10/20 14:07:59.55 hxLmqOB4.net
## 位数２ｎの二面体群の演算表を作る
# D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
n=7
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
e=1:n
d1=c(n,1:(n-1)) # 360°/n 回転移動
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i in 1:(n-1)){
Mnd[i+1,]=tikan(Mnd[i,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
172:卵の名無しさん
17/10/20 14:08:06.73 hxLmqOB4.net
t=n:1 ## 対称移動
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i in 1:(n-1)){
Mnt[i+1,]=tikan(Mnt[i,],d1)
}
rownames(Mnt)=c('t',paste0('td',1:(n-1)))
Mnt
(Mn=rbind(Mnd,Mnt))
names=rownames(Mn)
.L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1)
for(i in 1:(2*n)) .L[[i]]=Mn[i,]
names(.L)=names
.L
a=expand.grid(.L,.L)
(.a=mapply(g,a[,1],a[,2]))
(.M0=matrix(mapply(g,a[,1],a[,2]),length(.L)))
.M1=matrix(names[.a], length(.L))
rownames(.M1)=names
colnames(.M1)=names
.M1
print(.M1, quote=FALSE)
173:卵の名無しさん
17/10/20 14:08:34.26 hxLmqOB4.net
> print(.M1, quote=FALSE)
e d1 d2 d3 d4 d5 d6 t td1 td2 td3 td4 td5 td6
e e d1 d2 d3 d4 d5 d6 t td1 td2 td3 td4 td5 td6
d1 d1 d2 d3 d4 d5 d6 e td6 t td1 td2 td3 td4 td5
d2 d2 d3 d4 d5 d6 e d1 td5 td6 t td1 td2 td3 td4
d3 d3 d4 d5 d6 e d1 d2 td4 td5 td6 t td1 td2 td3
d4 d4 d5 d6 e d1 d2 d3 td3 td4 td5 td6 t td1 td2
d5 d5 d6 e d1 d2 d3 d4 td2 td3 td4 td5 td6 t td1
d6 d6 e d1 d2 d3 d4 d5 td1 td2 td3 td4 td5 td6 t
t t td1 td2 td3 td4 td5 td6 e d1 d2 d3 d4 d5 d6
td1 td1 td2 td3 td4 td5 td6 t d6 e d1 d2 d3 d4 d5
td2 td2 td3 td4 td5 td6 t td1 d5 d6 e d1 d2 d3 d4
td3 td3 td4 td5 td6 t td1 td2 d4 d5 d6 e d1 d2 d3
td4 td4 td5 td6 t td1 td2 td3 d3 d4 d5 d6 e d1 d2
td5 td5 td6 t td1 td2 td3 td4 d2 d3 d4 d5 d6 e d1
td6 td6 t td1 td2 td3 td4 td5 d1 d2 d3 d4 d5 d6 e
174:卵の名無しさん
17/10/20 14:21:16.31 hxLmqOB4.net
これが抜けていた。
f <- function(x){
re=NULL
for(i in 1:length(.L)){
if(equal(.L[[i]],x)) return(i) # i for names[i]
}
}
g <- function(x,y){
f(tikan2(x,y)) # first y, then x makes names[i]
}
175:卵の名無しさん
17/10/20 14:30:53.07 hxLmqOB4.net
## 位数２ｎの二面体群の演算表を作る
# D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
n=7
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
equal <- function(x,y) sum((x-y)^2)==0
176:卵の名無しさん
17/10/20 14:32:25.05 hxLmqOB4.net
f <- function(x){
re=NULL
for(i in 1:length(.L)){
if(equal(.L[[i]],x)) return(i) # i for names[i]
}
}
g <- function(x,y){
f(tikan2(x,y)) # first y, then x makes names[i]
}
e=1:n
d1=c(n,1:(n-1)) # 360°/n 回転移動
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i in 1:(n-1)){
Mnd[i+1,]=tikan(Mnd[i,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
177:卵の名無しさん
17/10/20 14:32:49.69 hxLmqOB4.net
t=n:1 ## 対称移動
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i in 1:(n-1)){
Mnt[i+1,]=tikan(Mnt[i,],d1)
}
rownames(Mnt)=c('t',paste0('td',1:(n-1)))
Mnt
(Mn=rbind(Mnd,Mnt))
names=rownames(Mn)
.L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1)
for(i in 1:(2*n)) .L[[i]]=Mn[i,]
names(.L)=names
.L
a=expand.grid(.L,.L)
(.a=mapply(g,a[,1],a[,2]))
(.M0=matrix(mapply(g,a[,1],a[,2]),length(.L)))
.M1=matrix(names[.a], length(.L))
rownames(.M1)=names
colnames(.M1)=names
.M1
print(.M1, quote=FALSE)
178:卵の名無しさん
17/10/20 14:41:11.07 hxLmqOB4.net
>170-172で動作確認
179:卵の名無しさん
17/10/21 07:30:34.81 H3fe+vIj.net
n=10での演算表
e d1 d2 d3 d4 d5 d6 d7 d8 d9 t td1 td2 td3 td4 td5 td6 td7 td8 td9
e e d1 d2 d3 d4 d5 d6 d7 d8 d9 t td1 td2 td3 td4 td5 td6 td7 td8 td9
d1 d1 d2 d3 d4 d5 d6 d7 d8 d9 e td9 t td1 td2 td3 td4 td5 td6 td7 td8
d2 d2 d3 d4 d5 d6 d7 d8 d9 e d1 td8 td9 t td1 td2 td3 td4 td5 td6 td7
d3 d3 d4 d5 d6 d7 d8 d9 e d1 d2 td7 td8 td9 t td1 td2 td3 td4 td5 td6
d4 d4 d5 d6 d7 d8 d9 e d1 d2 d3 td6 td7 td8 td9 t td1 td2 td3 td4 td5
d5 d5 d6 d7 d8 d9 e d1 d2 d3 d4 td5 td6 td7 td8 td9 t td1 td2 td3 td4
d6 d6 d7 d8 d9 e d1 d2 d3 d4 d5 td4 td5 td6 td7 td8 td9 t td1 td2 td3
d7 d7 d8 d9 e d1 d2 d3 d4 d5 d6 td3 td4 td5 td6 td7 td8 td9 t td1 td2
d8 d8 d9 e d1 d2 d3 d4 d5 d6 d7 td2 td3 td4 td5 td6 td7 td8 td9 t td1
d9 d9 e d1 d2 d3 d4 d5 d6 d7 d8 td1 td2 td3 td4 td5 td6 td7 td8 td9 t
t t td1 td2 td3 td4 td5 td6 td7 td8 td9 e d1 d2 d3 d4 d5 d6 d7 d8 d9
td1 td1 td2 td3 td4 td5 td6 td7 td8 td9 t d9 e d1 d2 d3 d4 d5 d6 d7 d8
td2 td2 td3 td4 td5 td6 td7 td8 td9 t td1 d8 d9 e d1 d2 d3 d4 d5 d6 d7
td3 td3 td4 td5 td6 td7 td8 td9 t td1 td2 d7 d8 d9 e d1 d2 d3 d4 d5 d6
td4 td4 td5 td6 td7 td8 td9 t td1 td2 td3 d6 d7 d8 d9 e d1 d2 d3 d4 d5
td5 td5 td6 td7 td8 td9 t td1 td2 td3 td4 d5 d6 d7 d8 d9 e d1 d2 d3 d4
td6 td6 td7 td8 td9 t td1 td2 td3 td4 td5 d4 d5 d6 d7 d8 d9 e d1 d2 d3
td7 td7 td8 td9 t td1 td2 td3 td4 td5 td6 d3 d4 d5 d6 d7 d8 d9 e d1 d2
td8 td8 td9 t td1 td2 td3 td4 td5 td6 td7 d2 d3 d4 d5 d6 d7 d8 d9 e d1
td9 td9 t td1 td2 td3 td4 td5 td6 td7 td8 d1 d2 d3 d4 d5 d6 d7 d8 d9 e
180:卵の名無しさん
17/10/21 19:00:23.36 TcIT7+c6.net
## D6で位数3の部分群を虱潰しに探す
a=t(combn(1:12,3))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
# names(H)=names(G[c(a[ii,1],a[ii,2],a[ii,3])])
M=matrix(NA,length(H),length(G))
for(i in 1:length(H)){
for(j in 1:length(G)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
#colnames(M)=paste0(names(G),'*')
#rownames(M)=names(H)
#print(M,quote=FALSE)
re=matrix(NA,nrow(M),ncol(M))
for(i1 in 1:ncol(M)){
re[,i1]=sort(M[,i1])
}
#print(re,quote=FALSE)
#print(unique(t(re)),quote=FALSE)
RES[ii] =nrow(unique(t(re)))
}
RES
idx=which(RES==4)
(b=a[idx,])
print(matrix(names[b],ncol=3),quote=FALSE)
181:卵の名無しさん
17/10/21 19:00:49.55 TcIT7+c6.net
> (b=a[idx,])
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
[3,] 7 9 11
[4,] 8 10 12
> print(matrix(names[b],ncol=3),quote=FALSE)
[,1] [,2] [,3]
[1,] e d2 d4
[2,] d1 d3 d5
[3,] t td2 td4
[4,] td1 td3 td5
182:卵の名無しさん
17/10/21 20:08:08.06 TcIT7+c6.net
## 二面体群Dnで位数mの部分群を虱潰しに探す
n=9 ; m=3
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f1 <- function(x,L){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(names(L)[i]) # i for names[i]
}
}
183:卵の名無しさん
17/10/21 20:19:32.79 TcIT7+c6.net
URLﾘﾝｸ(i.imgur.com)
184:卵の名無しさん
17/10/21 20:20:03.91 TcIT7+c6.net
##
a=t(combn(1:(2*n),m))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list() # H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
for(jj in 1:m){
H[[jj]]=G[[a[ii,jj]]]
}
M=matrix(NA,m,2*n)
for(i in 1:m){
for(j in 1:(2*n)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
re=matrix(NA,m,2*n)
for(i1 in 1:(2*n)){
re[,i1]=sort(M[,i1])
}
RES[ii] =nrow(unique(t(re)))
}
idx=which(RES==(2*n/m)) # Lagrangea's theorem
(b=a[idx,])
print(matrix(names[b],ncol=m),quote=FALSE)
185:卵の名無しさん
17/10/21 20:21:10.03 TcIT7+c6.net
> (b=a[idx,])
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
[4,] 10 13 16
[5,] 11 14 17
[6,] 12 15 18
> print(matrix(names[b],ncol=m),quote=FALSE)
[,1] [,2] [,3]
[1,] e d3 d6
[2,] d1 d4 d7
[3,] d2 d5 d8
[4,] t td3 td6
[5,] td1 td4 td7
[6,] td2 td5 td8
# 御明算
186:卵の名無しさん
17/10/21 21:11:33.21 TcIT7+c6.net
>>179
(debugged)
a=t(combn(1:(2*n),m))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list() # H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
for(jj in 1:m){
H[[jj]]=G[[a[ii,jj]]]
}
M=matrix(NA,m,2*n)
for(i in 1:m){
for(j in 1:(2*n)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
re=matrix(NA,m,2*n)
for(i1 in 1:(2*n)){
re[,i1]=sort(M[,i1])
}
RES[ii] =nrow(unique(t(re)))
}
idx=which(RES==(2*n/m)) # Lagrangea's theorem
b=a[idx,]
(B=b[b[,1]==1,])
print(matrix(names[B],ncol=m),quote=FALSE)
187:卵の名無しさん
17/10/21 23:55:34.83 TcIT7+c6.net
## 二面体群Dnですべての部分群を返す
dihedral_group <- function(n=9,m=3){ # start
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f1 <- function(x,L){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(names(L)[i]) # i for names[i]
}
}
g12 <- function(x,y,L=G){
f1(tikan2(x,y),L)
}
equal <- function(x,y) sum((x-y)^2)==0
188:卵の名無しさん
17/10/21 23:56:01.41 TcIT7+c6.net
e=1:n
d1=c(n,1:(n-1)) # rotation
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i0 in 1:(n-1)){
Mnd[i0+1,]=tikan(Mnd[i0,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
t=n:1 ## mirror
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i2 in 1:(n-1)){
Mnt[i2+1,]=tikan(Mnt[i2,],d1)
}
189:Mnt rownames(Mnt)=c('t',paste0('td',1:(n-1))) Mn=rbind(Mnd,Mnt) names=rownames(Mn) Mn .L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1) for(i3 in 1:(2*n)) .L[[i3]]=Mn[i3,] names(.L)=names G=.L
190:卵の名無しさん
17/10/21 23:56:34.66 TcIT7+c6.net
a=t(combn(1:(2*n),m))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list() # H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
for(jj in 1:m){
H[[jj]]=G[[a[ii,jj]]]
}
M=matrix(NA,m,2*n)
for(i in 1:m){
for(j in 1:(2*n)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
re=matrix(NA,m,2*n)
for(i1 in 1:(2*n)){
re[,i1]=sort(M[,i1])
}
RES[ii] =nrow(unique(t(re)))
}
RES
idx=which(RES==(2*n/m)) # Lagrangea's theorem
(b=a[idx,])
#b=as.matrix(b)
if(m==1||m==2*n){
B=1:m
}else{
B=b[b[,1]==1,]
}
matrix(names[B],ncol=m)
} # End of Function, dihedral_group
191:卵の名無しさん
17/10/21 23:57:03.96 TcIT7+c6.net
.print <- function(x) print(x,quote=FALSE)
sub_group <- function(n){
N=2*n
xx=which(N%%(1:N)==0)
res=lapply(xx,function(x)dihedral_group(n,x))
names(res)=xx
.print(res)
}
192:卵の名無しさん
17/10/21 23:57:45.77 TcIT7+c6.net
> sub_group(6)
$`1`
[,1]
[1,] e
$`2`
[,1] [,2]
[1,] e d3
[2,] e t
[3,] e td1
[4,] e td2
[5,] e td3
[6,] e td4
[7,] e td5
$`3`
[,1] [,2] [,3]
[1,] e d2 d4
$`4`
[,1] [,2] [,3] [,4]
[1,] e d3 t td3
[2,] e d3 td1 td4
[3,] e d3 td2 td5
$`6`
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] e d1 d2 d3 d4 d5
[2,] e d2 d4 t td2 td4
[3,] e d2 d4 td1 td3 td5
$`12`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] e d1 d2 d3 d4 d5 t td1 td2 td3 td4 td5
193:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 09:18:14.91 bCTtjSj1.net
# Greatest Common Divisor with Euclidean algorithm
GCD <- function(x,y){
if(round(x)!=x | round(y)!=y) stop('Not integer')
if(!x&!y) return(NA)
a=max(c(abs(x),abs(y)))
b=min(c(abs(x),abs(y)))
if(!a*b) return(a)
r=integer()
r[1]=a
r[2]=b
i=1
r[i+2]=r[i]%%r[i+1]
while(r[i+2]!=0 & r[i+2]!=1){
i=i+1
r[i+2]=r[i]%%r[i+1]
}
return(ifelse(r[i+2]==0,r[i+1],1))
}
GCD(2.3,4)
GCD(0,0)
GCD(24,15)
GCD(16,15)
GCD(5,0)
GCD(24,-15)
GCD(-24,-15)
194:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 13:48:21.69 bCTtjSj1.net
> sub_group(9)
$`1`
[,1]
[1,] e
$`2`
[,1] [,2]
[1,] e t
[2,] e td1
[3,] e td2
[4,] e td3
[5,] e td4
[6,] e td5
[7,] e td6
[8,] e td7
[9,] e td8
$`3`
[,1] [,2] [,3]
[1,] e d3 d6
$`6`
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] e d3 d6 t td3 td6
[2,] e d3 d6 td1 td4 td7
[3,] e d3 d6 td2 td5 td8
$`9`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] e d1 d2 d3 d4 d5 d6 d7 d8
$`18`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] e d1 d2 d3 d4 d5 d6 d7 d8 t td1 td2 td3 td4 td5
[,16] [,17] [,18]
[1,] td6 td7 td8
195:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 15:51:11.54 bCTtjSj1.net
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
# Gを有限群とすると g∈G n=|G| （位数） g^n=e :単位元
n=length(G)
re=list()
for(i in 1:n){
tikan2(G[[i]],G[[i]])
for(j in 1:n) data[[j]]=G[[i]]
re[[i]]=Reduce(tikan2,data) # Reduce(function(x,y)x+y,c(1,2,3,4),accumulate = TRUE)
}
> re
[[1]]
[1] 1 2 3 4 5 6 7 8 9
[[2]]
[1] 1 2 3 4 5 6 7 8 9
......
[[17]]
[1] 1 2 3 4 5 6 7 8 9
[[18]]
[1] 1 2 3 4 5 6 7 8 9
196:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 16:46:17.03 bCTtjSj1.net
> Reduce(function(x,y)x+y,c(1,2,3,4,5),accumulate = TRUE)
[1] 1 3 6 10 15
> Reduce('+',c(1,2,3,4,5),accu=TRUE) ; cumsum(1:5)
[1] 1 3 6 10 15
[1] 1 3 6 10 15
> Reduce('*',c(1,2,3,4,5)) ; factorial(5)
[1] 120
[1] 120
> beki2 <- function(x,y,p,...){ #x^y%%p
+ if(y==0) return(1%%p)
+ data=rep(x,y)
+ Reduce(function(x,y) (x*y)%%p, data,...)
+ }
> beki2(2,10,10,accumulate=TRUE)
[1] 2 4 8 6 2 4 8 6 2 4
> beki2(2,100,10)
[1] 6
> beki2(2,0,10)
[1] 1
197:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 17:29:41.55 bCTtjSj1.net
# When p is a prime number and a is coprime to p,
#
# a^(p-1)≡1 (mod p)
#
# Check if it works when p is lower than N.
N=100
(p100=(1:N)[!!gmp::isprime(1:N)] ) # primes < 100
m=length(p100)
cop=list()
for(i in 1:m){
p=p100[i]
cop[[i]]=(1:p)[gmp::gcd(1:p,p)==1] # disjoint,coprime to p100[i]
}
names(cop)=p100
cop
re=cop
for(i in 1:m){
p=p100[i]
for(j in 1:length(cop[[i]])){
re[[i]][j] = beki2(cop[[i]][j],p-1,p) # coprime^(p-1)%%p
}
}
re
> str(re)
List of 25
$ 2 : int 1
$ 3 : int [1:2] 1 1
...
$ 89: int [1:88] 1 1 1 1 1 1 1 1 1 1 ...
$ 97: int [1:96] 1 1 1 1 1 1 1 1 1 1 ...
フェルマーの定理
198:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 18:21:14.08 r/LJPTCR.net
>>191
フェルマーも�
199:ｬ定理の方ね
200:卵の名無しさん
17/10/25 19:11:29.66 nTDNW+TM.net
>>122
##
n=10
conf.level=0.95
PMF <- function(x) x^n # AUC==integrate(PMF,0,1) == 1/(n+1)
PDF <- function(x) (n+1)*x^n # PDF==PMF/AUC
Ex= (n+1)/(n+2) # integrate(x*PDF,0,1)==integrate((n+1)*x^(n+1),0,1)== (n+1)/(n+2)
CDF <- function(x) x^(n+1) # == integrate(PDF,0,x)
lwr = (1-conf.level)^(1/(n+1)) # CDF(lwr) == 1-conf.level
exp(log(1-conf.level)/(n+1))
curve(PDF)
curve(x*PDF(x))
integrate(function(x) x*PDF(x),0,1)$value
integrate(function(x) x*(n+1)*x^n,0,1)$value # integrate(function(x)(n+1)*x^(n+1),0,1)
(n+1)/(n+2)
# |((n+1)/(n+2))*x^(n+2)|[0,1]
201:卵の名無しさん
17/10/25 19:13:22.53 nTDNW+TM.net
>>193
## binomとの比較
binom::binom.confint(c(1,10,100),c(1,10,100),conf=0.95)
# n out of n
# mean=(n+1)/(n+2)
# lower=(n+1)√(1-conf.level) = 0.05^(1/(n+1))
non=function(n,conf.level=0.95){ # n hits out of n shoots
b=binom::binom.confint(n,n)[,c('method','mean','lower')]
n_1_root=function(n)(0.05)^(1/(n+1))
a=data.frame(method='root(n+1)',mean=(n+1)/(n+2),lower=(1-conf.level)^(1/(n+1)))
rbind(a,b)
}
> non(1)
method mean lower
1 root(n+1) 0.6666667 0.22360680
2 agresti-coull 1.0000000 0.16749949
3 asymptotic 1.0000000 1.00000000
4 bayes 0.7500000 0.22851981
5 cloglog 1.0000000 0.02500000
6 exact 1.0000000 0.02500000
7 logit 1.0000000 0.02500000
8 probit 1.0000000 0.02500000
9 profile 1.0000000 0.13811125
10 lrt 1.0000000 0.14650325
11 prop.test 1.0000000 0.05462076
12 wilson 1.0000000 0.20654931
202:卵の名無しさん
17/10/25 19:13:48.14 nTDNW+TM.net
> non(10)
method mean lower
1 root(n+1) 0.9166667 0.7615958
2 agresti-coull 1.0000000 0.6791127
3 asymptotic 1.0000000 1.0000000
4 bayes 0.9545455 0.8292269
5 cloglog 1.0000000 0.6915029
6 exact 1.0000000 0.6915029
7 logit 1.0000000 0.6915029
8 probit 1.0000000 0.6915029
9 profile 1.0000000 0.7303058
10 lrt 1.0000000 0.8252466
11 prop.test 1.0000000 0.6554628
12 wilson 1.0000000 0.7224672
> non(100)
method mean lower
1 root(n+1) 0.9901961 0.9707748
2 agresti-coull 1.0000000 0.9555879
3 asymptotic 1.0000000 1.0000000
4 bayes 0.9950495 0.9810231
5 cloglog 1.0000000 0.9637833
6 exact 1.0000000 0.9637833
7 logit 1.0000000 0.9637833
8 probit 1.0000000 0.9637833
9 profile 1.0000000 0.9670434
10 lrt 1.0000000 0.9809757
11 prop.test 1.0000000 0.9538987
12 wilson 1.0000000 0.9630065
203:卵の名無しさん
17/10/25 19:24:32.39 nTDNW+TM.net
>>195
n=10 の時のPDFと９５％CI下限値0.7615958
URLﾘﾝｸ(i.imgur.com)
204:卵の名無しさん
17/10/25 19:24:57.43 nTDNW+TM.net
>>195
n=10 の時のPDFと９５％CI下限値0.7615958
URLﾘﾝｸ(i.imgur.com)
205:卵の名無しさん
17/10/25 20:40:32.01 nTDNW+TM.net
#.n発r中の狙撃手が.N発狙撃するときの命中数を返す
Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
print(summary(SS))
invisible(SS)
}
206:卵の名無しさん
17/10/26 12:27:43.35 ljbBi1cA.net
Golgo <- function(n=3,r=1,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
AUC=integrate(PMF,0,1)$value
# library(hypergeo)
# f.auc <- function(x) choose(n,r)*x^(r+1)/(r+1)*hypergeo(r+1,r-n,r+2,x)
# auc=as.numeric(f.auｎc(1)-f.auc(0)) ; 1/auc
PDF <- function(x)PMF(x)/AUC
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),1),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
207:卵の名無しさん
17/10/26 16:21:26.13 ljbBi1cA.net
# ｆ回失敗した後にｓ回目の成功,K:シミュレーション回数
Dotsubo <- function(s=1,f=4,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(p) choose(s+f-1,f)*(1-p)^f*p^s
AUC= (gamma(s+f)*gamma(s+1)) / (gamma(s)*gamma(f+s+2))
PDF <- function(p) (1-p)^f*p^s*gamma(f+s+2)/(gamma(f+1)*gamma(s+1))
curve(PDF)
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),2),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
208:卵の名無しさん
17/10/26 17:43:49.64 z4lFiDnk.net
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
209:卵の名無しさん
17/10/27 15:54:12.90 dzxKDqmi.net
# ｎ発ｒ中の狙撃手がＮ発狙撃するときの命中数を返す
Golgo.sim <- function(N, n, r, k=10^3,Print=TRUE){ # k:シミュレーション回数
f <-function(S,.N=N,.n=n){ # 成績サンプル：命中S個、外れ(N-S)個
y=c(rep(1,S),rep(0,.N-S))
sum(sample(y,.n)) # その成績サンプルからｎ個数取り出したときの命中数
}
xx=r:N # ｒ未満ではｒ個命中することはないので除外
SS=NULL # 容れ子
for(i in 1:k){
x=sapply(xx,f) # 命中数の配列
SS=append(SS,which(x==r)-1+r) # 命中数がｒ個のときの成績サンプルの命中数Ｓの配列をつくる
}
print(summary(SS))
print(quantile(SS,probs=c(0.025,0.05,0.95,0.975)))
print(c(mode=names(which.max(table(SS)))),quote=FALSE)
if(Print) {
hist(SS,xlim=c(0,N),freq=FALSE,col=sample(colors(),1),main='',xlab='Hits')
lines(density(SS))}
invisible(SS)
}
210:卵の名無しさん
17/10/27 15:55:45.80 pwjeiI6l.net
# n発r中の期待値
Golgo <- function(n=3,r=1,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
AUC=integrate(PMF,0,1)$value
# library(hypergeo)
# f.auc <- function(x) choose(n,r)*x^(r+1)/(r+1)*hypergeo(r+1,r-n,r+2,x)
# auc=as.numeric(f.auc(1)-f.auc(0))
PDF <- function(x)PMF(x)/AUC
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),1),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
211:卵の名無しさん
17/10/29 17:45:50.97 LlbU36d2.net
data{// coin5.stan
int N;
int<lower=0,upper=1> Y[N];
}
parameters{
real<lower=0,upper=1> p;
}
model{
for(n in 1:N)
Y[n] ~ bernoulli(p);
}
data{// coin1.stan
int<lower=0,upper=1> Y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
Y ~ bernoulli(p);
}
212:卵の名無しさん
17/10/29 22:16:10.95 LlbU36d2.net
# クラインの4元群Vが正6面体群S(P6)の正規部分群であることの確認
equal <- function(x,y) sum((x-y)^2)==0
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f2 <- function(x,L=G){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(i) # i for names[i]
}
}
a =c(2,1,4,3)
b =c(3,4,1,2)
c =c(4,3,2,1)
e =c(1,2,3,4)
d =c(1,3,4,2)
t =c(1,2,4,3)
d2 =tikan2(d,d) # 1 4 2 3 # d3=tikan2(d2,d) == e
td =tikan2(t,d) # 1 3 2 4
td2=tikan2(td,d)
213:卵の名無しさん
17/10/29 22:17:23.29 LlbU36d2.net
V=list(e,a,b,c)
D3=list(e,d,d2,t,td,td2)
gr=expand.grid(V,D3)
.m=mapply(tikan2,gr[,1],gr[,2])
t.m=t(.m)
G=list()
for(i in 1:nrow(t.m)) G[[i]]=t.m[i,]
names(G)=paste0('t',1:length(G))
G # G: S(P6)
lG=length(G)
V
lV=length(V)
.M1=matrix(NA,nrow=lV,ncol=lG)
for (i in 1:lG){
for(j in 1:lV){
.M1[j,i]=f2(tikan2(G[[i]],V[[j]]))
}
}
.M2=matrix(NA,nrow=lV,ncol=lG)
for (i in 1:lG){
for(j in 1:lV){
.M2[j,i]=f2(tikan2(V[[j]],G[[i]]))
}
}
identical(apply(.M1,2,sort) ,apply(.M2,2,sort))
> identical(apply(.M1,2,sort) ,apply(.M2,2,sort))
[1] TRUE
214:卵の名無しさん
17/10/30 07:00:24.65 OjMfeTM6.net
## 写像を元とする集合の積を計算する
equal <- function(x,y) sum((x-y)^2)==0 & length(x)==length(y)
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
witch <- function(x,L=G){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(i) # i for names[i]
}
}
tikan3 <- function(X,Y,L=G){ # replace Y first then X and return index of L
Z=tikan2(X,Y)
witch(Z,L)
}
product <- function(A,B,L=G){ # product of sets, retur index of L
la=length(A)
lb=length(B
215:) gr=expand.grid(1:la,1:lb) f<-function(x,y) tikan3(A[[x]],B[[y]],L) re=mapply(f,gr[,1],gr[,2]) return(sort(unique(re))) }
216:卵の名無しさん
17/10/30 07:01:59.31 OjMfeTM6.net
GROUP <- function(n=6){
e=1:n
d1=c(n,1:(n-1)) # rotation
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i0 in 1:(n-1)){
Mnd[i0+1,]=tikan(Mnd[i0,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
t=n:1 ## mirror
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i2 in 1:(n-1)){
Mnt[i2+1,]=tikan(Mnt[i2,],d1)
}
Mnt
rownames(Mnt)=c('t',paste0('td',1:(n-1)))
Mn=rbind(Mnd,Mnt)
names=rownames(Mn)
Mn
.L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1)
for(i3 in 1:(2*n)) .L[[i3]]=Mn[i3,]
names(.L)=names
G=.L
return(G)
}
217:卵の名無しさん
17/10/31 16:10:43.91 1m4iMpv8.net
HDCI <- function(PMF,cl=0.95){ # Highest Density Confidence Interval
PDF=PMF/sum(PMF)
rsPDF=rev(sort(PDF))
min.density=rsPDF[min(which(cumsum(rsPDF)>=cl))]
index=which(PDF>=min.density)
data.frame(lower.idx=round(min(index)),upper.idx=round(max(index)),actual.CI=sum(PDF[index]))
}
218:卵の名無しさん
17/10/31 16:11:33.50 1m4iMpv8.net
# n発r中の期待値
Golgo2 <- function(n=3,r=1,cl=0.95,k=0.00001,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
xx=seq(0,1,by=k)
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
print(c(lower=lower,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=5)
if(Print){plot(xx,pmf,lwd=1,xlab='probability of hit',
ylab=paste('probabilty of',r,'hits out of ',n,'shoots'))}
}
219:卵の名無しさん
17/10/31 20:38:26.69 A7oOP/mA.net
# １年：進級失敗１０人、うち１人放校
# ２年：進級失敗１６人、放校なし
# ３年：進級失敗３４人、うち放校９人
# ４年：進級失敗９人、うち放校２人
# ５年：進級失敗１０人、うち放校１人
# ６年：卒業失敗２６人、うち放校１人
# 一学年約120～130人前後。
#ｽﾚﾘﾝｸ(doctor板:1番)
flunk=c(10,16,34,9,10,26)
expel=c(1,0,9,2,1,1)
n=125
total.fl=sum(flunk)
total.ex=sum(expel)
re=matrix(NA,6,1)
for(i in 1:6){
re[i,]=poisson.test(c(flunk[i],n) ,c(total.fl,n*6))$p.value
}
re
for(i in 1:6){
re[i,]=poisson.test(c(expel[i],n) ,c(total.ex,n*6))$p.value
}
re
prop.test(flunk,rep(n,6))
pairwise.prop.test(flunk,rep(n,6),'bon')
prop.test(expel,rep(n,6))
d=cbind(expel,rep(n,6)-expel) ; d
fisher.test(d)
pairwise.prop.test(expel,rep(n,6),'none')
pairwise.prop.test(expel,rep(n,6),'bon')
220:卵の名無しさん
17/11/01 09:28:50.55 zVXhNw2e.net
## 確率分布から信頼区間を出す
HDCI2 <- function(PMF,cl=0.95,k=0.0001,Print=TRUE){
xx=seq(0,1,by=k)
xx=xx[-1]
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
print(c(lower=lower,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=4)
if(Print) plot(xx,pmf,xlab='prior',ylab='posterior')
}
HDCI2(function(x)dnbinom(5-1,1,x))
221:卵の名無しさん
17/11/01 09:38:36.29 zVXhNw2e.net
## 確率分布から最尤値・期待値・信頼区間を出す
HDCI2 <- function(PMF,cl=0.95,k=0.0001,Print=FALSE){
pp=seq(0,1,by=k)
xx=pp[-1]
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
mode=xx[which.max(pmf)]
print(c(lower=lower,mode=mode,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=4)
if(Print) plot(xx,pmf,xlab='prior',ylab='posterior')
}
222:卵の名無しさん
17/11/05 17:08:17.51 G4ZpDCNG.net
NHST Has 100% False Alarm Rate in Sequential
Testing
Under NHST, sequential testing of data generated from the null
hypothesis will eventually lead to a false alarm. With infinite
patience, there is 100% probability of falsely rejecting the null.
This is known as “sampling to reach a foregone conclusion” (e.g.,
Anscombe, 1954). To illustrate this phenomenon, a computer
simulation generated random values from a normal distribution
with mean zero and standard deviation one, assigning each sequential
value alternately to one or the other of two groups, and at each
step conducting a two-group t test assuming the current sample
sizes were fixed in advance. Each simulat
223:ed sequence began with N1 N2 3. If at any step the t test indicated p .05, the sequence was stopped and the total N ( N1 N2) was recorded. Otherwise, another random value was sampled from the zero-mean normal and included in the smaller group, and a new t test was conducted. For purposes of illustration, each sequence was limited to a maximum sample size of N 5,000. The simulation ran 1,000 sequences. NHST <- function(N=5000){ x=rnorm(3) y=rnorm(3) while(length(x) < N & t.test(x,y,var.equal=TRUE)$p.value >= 0.05){ x=append(x,rnorm(1)) y=append(y,rnorm(1)) } return(length(x)) } res <- replicate(1000,NHST())
224:卵の名無しさん
17/11/05 17:27:50.94 G4ZpDCNG.net
# NHST Has 100% False Alarm Rate in Sequential Testing (NHST : Null Hypothesis Significance Testing)
# Under NHST, sequential testing of data generated from the null
# hypothesis will eventually lead to a false alarm. With infinite
# patience, there is 100% probability of falsely rejecting the null.
# This is known as “sampling to reach a foregone conclusion” (e.g.,
# Anscombe, 1954). To illustrate this phenomenon, a computer
# a computer simulation generated random values from a normal distribution
# with mean zero and standard deviation one, assigning each sequential
# value alternately to one or the other of two groups, and at each
# step conducting a two-group t test assuming the current sample
# sizes were fixed in advance. Each simulated sequence began with
# N1=N2=3. If at any step the t test indicated p .05, the
# sequence was stopped and the total N (=N1+N2) was recorded.
NHST <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & t.test(x,y,var.equal=TRUE)$p.value >= 0.05){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
FRR <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
NN=c(10,20,40,80,160,320,640,1280,2560,5120,10240)
RjR=sapply(NN,FRR)
plot(NN,RjR,type='b')
225:卵の名無しさん
17/11/05 18:24:17.31 G4ZpDCNG.net
# Bayesian decision making, using the HDI and ROPE, does not
# suffer a 100% false alarm rate in sequential testing. Instead, the
# false alarm rate asymptotes at a much lower level, depending on
# the choice of ROPE. For illustration, again a computer simulation
# generated random values from a normal distribution with mean of
# zero and standard deviation of one, assigning each sequential value
# alternately to one or the other of two groups but at each step
# conducting a Bayesian analysis and checking whether the 95%
# HDI completely excluded or was contained within a ROPE from 0.15 to 0.15
# ROPE* region of practical equivalence
librry(BEST)
BA <- function(N){
mc=BEST::BESTmcmc(rnorm(N),rnorm(N),numSavedSteps=1000, burnInSteps=500)
re=summary(mc,ROPEm=c(-0.15,15))
return(re[['muDiff','%InROPE']])
}
NN=c(10,20,40,80,160,320,640,1280,2560,5120)
inrope=sapply(NN,BA)
plot(NN,inrope,type='b',pch=19,xlab='N',ylab='% in ROPE')
226:卵の名無しさん
17/11/05 20:18:58.81 G4ZpDCNG.net
Bayesian Estimation Supersedes the t-test (BEST) - online
URLﾘﾝｸ(www.sumsar.net)
トレースプロットがみられるのがうれしい。
227:卵の名無しさん
17/11/05 20:22:56.18 G4ZpDCNG.net
新版の予約受付中か。
URLﾘﾝｸ(www.amazon.co.jp)
ベイズ統計がどれくらい組み込まれたをみてから買うかな。
228:卵の名無しさん
17/11/05 20:28:14.21 G4ZpDCNG.net
# A君の彼女は女子大生、B君の彼女は女子高生。
# Y1女子大生n1=100人とY2女子高生n2=100人の胸囲を測定して
# 前者が平均82 , 標準偏差3
# 後者が平均81 , 標準偏差3
# 有意差はあるか？
T.test=function(n1,n2,m1,m2,sd1,sd2){
SE12=sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
T=(m1-m2)/SE12
2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
T.test(100,100,82,81,3,3)
y1=82+scale(rnorm(100))*3
y2=81+scale(rnorm(100))*3
t.test(y1,y2,var.equal = TRUE)
library(BEST)
BESTout <- BESTmcmc(y1,y2)
plot(BESTout,ROPE=c(-2,2))
summary(BESTout,ROPEm=c(-2,2))
plotAll(BESTout,ROPEm=c(-2,2))
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-0.5,0.5), ROPEsd=c(-1,1), ROPEeff=c(-0.5,0.5),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=5)
powerPro
229:卵の名無しさん
17/11/05 20:33:01.04 G4ZpDCNG.net
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-2,2), ROPEsd=c(-0,0), ROPEeff=c(-0,0),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=1000)
powerPro
230:卵の名無しさん
17/11/05 20:51:07.47 G4ZpDCNG.net
Jikkan <- function(x){
re=summary(BESTout,ROPEm=c(-x,x))
re[['muDiff','%InROPE']]
}
xx=seq(0,2.5,by=0.01)
InRope=sapply(xx,Jikkan)
plot(xx,InRope,type='l',lwd=2, xlab='Difference',ylab='% Practically Equivalent')
abline(h=95, lty=3)
(Diff=uniroot(function(x,u0=95) Jikkan(x)-u0, c(1.5,2))$root)
plot(BESTout,ROPE=c(-Diff,Diff))
URLﾘﾝｸ(i.imgur.com)
231:卵の名無しさん
17/11/06 21:07:39.66 1uZMaFLW.net
##
library(BayesFactor)
tBF <- function(x,y){
t=t.test(x,y,var.equal = TRUE)$statistic
ttest.tstat(t,length(x),length(y),simple=TRUE)
}
NHST2 <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & tBF(x,y) <= 1){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
FRR2 <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST2(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
> FRR2(10)
[1] 0.37
> FRR2(50)
[1] 0.46
> FRR2(100)
[1] 0.56
ベイズ因子をつかっても ｔ 検定を組み込むと
# NHST Has 100% False Alarm Rate in Sequential Testing (NHST : Null Hypothesis Significance Testing)
という結論は変わらないな。
232:卵の名無しさん
17/11/06 21:51:42.88 1uZMaFLW.net
>>216
Sequential Samplingに修正
library(BEST)
BA <- function(x,y){
mc=BEST::BESTmcmc(x,y,num=5000)
re=summary(mc,ROPEm=c(-0.1,0.1))
return(re[['muDiff','%InROPE']])
}
x=rnorm(3)
y=rnorm(3)
inROPE=numeric()
N=100
for(i in 1:(N-3)){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
inROPE[i]=BA(x,y)
}
plot(3:N,inROPE,type='b',pch=19,xlab='n',ylab='% in ROPE')
233:卵の名無しさん
17/11/07 00:05:34.55 2r/5Mjhu.net
>>201
n=1;r=1
f=function(p)choose(n,r)*p^r*(1-p)^(n-r)
AUC=integrate(f,0,1)$value
integrate(function(x) x*f(x)/AUC,0,1)$value
234:卵の名無しさん
17/11/07 00:49:20.56 2r/5Mjhu.net
library(BayesFactor)
tBF2 <- function(x,y,...){
tt=t.test(x,y,var.equal = TRUE)
bf=ttest.tstat(tt$statistic,length(x),length(y),simple=TRUE,...)
p=tt$p.value
return(c(bf=1/bf,p=p))
}
N=10
BF=numeric()
p.value=numeric()
k=1000
for(i in 1:k){
x=rnorm(N) ; y=rnorm(N)
bfp=tBF2(x,y)
BF[i]=bfp[1]
p.value[i]=bfp[2]
}
.m=cbind(BF,p.value)
.m[which(p.value<0.05),]
mean(p.value<0.05)
.m[which(BF<1),]
mean(BF<1)
235:卵の名無しさん
17/11/07 20:40:25.24 2r/5Mjhu.net
# crooked coin or dice
crooked <- function(n,r,H0=0.5,d=0.1,cl=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,cl,shape1=1+r,shape2=1+n-r)
if(Print){
ROPEl=H0*(1-d)
ROPEu=H0*(1+d)
curve(dbeta(x,1+r,1+n-r),lwd=1,xlab='p',ylab='density',...) # posterior
segments(ROPEl,0,ROPEu,lwd=5,col='navy')
segments(hdi[1],0,hdi[1],dbeta(hdi[1],1+r,1+n-r),col='blue')
segments(hdi[2],0,hdi[2],dbeta(hdi[2],1+r,1+n-r),col='blue')
legend('topleft',bty='n',legend=c('Range of Practival Equivalence',
'High Density Interval'),lwd=c(5,1),col=c('navy','blue'))
}
return(hdi)
}
> crooked(3000,490,H0=1/6,xlim=c(0.1,0.3))
lower upper
0.1503971 0.1768440
236:卵の名無しさん
17/11/07 21:15:32.89 2r/5Mjhu.net
# crooked coin or dice
crooked <- function(n,r,H0=0.5,d=0.1,credMass=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,credMass,shape1=1+r,shape2=1+n-r)
if(Print){
ROPEl=H0*(1-d)
ROPEu=H0*(1+d)
curve(dbeta(x,1+r,1+n-r),lwd=1,xlab='p',ylab='density',...) # posterior
segments(ROPEl,0,ROPEu,lwd=5,col='navy')
segments(hdi[1],0,hdi[1],dbeta(hdi[1],1+r,1+n-r),col='blue')
segments(hdi[2],0,hdi[2],dbeta(hdi[2],1+r,1+n-r),col='blue')
legend('topleft',bty='n',legend=c('Range of Practival Equivalence',
paste0(credMass/0.01,'% High Density Interval')),lwd=c(5,1),col=c('navy','blue'))
}
return(hdi)
}
crooked(3000,490,H0=1/6,xlim=c(0.12,0.20))
237:卵の名無しさん
17/11/08 19:29:03.43 vjtLqicz.net
data{// coin5d.stan
int N;
int<lower=0,upper=1> Y[N];
real b; // border
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real ura = step(b-p); // ifelse(b-p>0,1,0)
}
model{
for(n in 1:N)
Y[n] ~ bernoulli(p);
}
238:卵の名無しさん
17/11/08 19:36:21.87 vjtLqicz.net
N=10
r=5
Y=c(rep(1,r),rep(0,N-r))
b=0.6
data <- list(Y=Y, N=N,b=b)
model.coin5d <- stan_model('coin5d.stan')
fit.coin5 <- sampling(model.coin5d, data=data, seed=123)
print(fit.coin5,digits=4)
Inference for Stan model: coin5d.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
p 0.4977 0.0035 0.1410 0.2286 0.3988 0.4966 0.5976 0.7744 1620 1.0036
ura 0.7552 0.0098 0.4300 0.0000 1.0000 1.0000 1.0000 1.0000 1933 1.0007
lp__ -8.8598 0.0234 0.7805 -11.0878 -9.0348 -8.5596 -8.3753 -8.3186 1116 1.0027
239:卵の名無しさん
17/11/08 20:57:34.40 vjtLqicz.net
# NHST Has 100% False Alarm Rate in Sequential Testing
# NHST : Null Hypothesis Significance Testing
# Under NHST, sequential testing of data generated from the null
# hypothesis will eventually lead to a false alarm. With infinite
# patience, there is 100% probability of falsely rejecting the null.
# This is known as “sampling to reach a foregone conclusion” (e.g.,
# Anscombe, 1954). To illustrate this phenomenon,
# a computer simulation generated random values from a normal distribution
# with mean zero and standard deviation one, assigning each sequential
# value alternately to one or the other of two groups, and at each
# step conducting a two-group t test assuming the current sample
# sizes were fixed in advance. Each simulated sequence began with
# N1=N2=3. If at any step the t test indicated p < .05, the
# sequence was stopped and the total N (=N1+N2) was recorded.
240:卵の名無しさん
17/11/08 20:57:42.57 vjtLqicz.net
FUN = wilcox.test
FUN = perm::permTS
FUN = lawstat::brunner.munzel.test
NHST <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & FUN(x,y)$p.value >= 0.05){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
NHST(100)
FRR <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
FRR(10)
FRR(50)
FRR(100)
FRR(500)
n=10 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
n=50 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
n=100 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
n=500 ; mean(replicate(k,FUN(rnorm(n),rnorm(n))$p.value) < 0.05)
241:卵の名無しさん
17/11/11 09:46:29.46 EghxJpRr.net
source('URLﾘﾝｸ(github.com)')
source('URLﾘﾝｸ(github.com)')
a=1;b=1
N=21
z=21
pp=seq(0,1,by=0.001)
likelihood=pp^z*(1-pp)^(N-z)
prior=dbeta(pp,a,b)
posterior = BernGrid(pp,prior/sum(prior),c(rep(1,z),rep(0,N-z)), plotType="Bars" ,
showCentTend="Mean" , showHDI=TRUE , showpD=TRUE )
242:卵の名無しさん
17/11/11 10:09:38.03 EghxJpRr.net
# core concept
par(mfrow=c(2,2))
z=1 ; N=1
pp=seq(0,1,by=0.001)
prior=pmin(pp,1-pp) # a isosceles triangular distribution
plot(pp,prior,type='h',lwd=5,col='lightblue')
likelihood=pp^z*(1-pp)^(N-z)
plot(pp,likelihood,type='h',col='lightblue')
(pD=sum(prior*likelihood))
posterior=prior*likelihood/pD
plot(pp,posterior,type='h', col=' maroon')
243:卵の名無しさん
17/11/11 11:35:44.48 ARZ4w8db.net
require(rjags)
z=21
N=21
y=c(rep(1,z),rep(0,N-z))
data=list(N=N,y=y)
modelString = "
model {
for (i in 1:N) {
y[i] ~ dbern(theta)
}
theta ~ dbeta(1,1)
}
"
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=data , # inits=initsList ,
n.chains=3 , n.adapt=500 )
update( jagsModel , n.iter=500 )
codaSamples = coda.samples( jagsModel , variable.names=c("theta") ,
n.iter=3334 )
summary(codaSamples)
plot(codaSamples)
244:卵の名無しさん
17/11/11 16:50:37.66 AQmX3Wf1.net
genMCMC <- # prior : theta ~ Beta(a,b)
function( data , numSavedSteps=50000 , saveName=NULL, a=1,b=1) {
require(rjags)
y = data$y
s = as.numeric(data$s) # converts character to consecutive integer levels
# Do some checking that data make sense:
if ( any( y!=0 & y!=1 ) ) { stop("All y values must be 0 or 1.") }
Ntotal = length(y)
Nsubj = length(unique(s))
# Specify the data in a list, for later shipment to JAGS:
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
#-----------------------------------------------------------------------------
# THE MODEL.
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern(
245: theta[s[i]] ) } for ( sIdx in 1:Nsubj ) { theta[sIdx] ~ dbeta(", a,',' , b," ) } } ") # close quote for modelString writeLines( modelString , con="TEMPmodel.txt" )
246:卵の名無しさん
17/11/11 16:50:57.13 AQmX3Wf1.net
#-----------------------------------------------------------------------------
# INTIALIZE THE CHAINS.
# Initial values of MCMC chains based on data:
# Option 1: Use single initial value for all chains:
# thetaInit = rep(0,Nsubj)
# for ( sIdx in 1:Nsubj ) { # for each subject
# includeRows = ( s == sIdx ) # identify rows of this subject
# yThisSubj = y[includeRows] # extract data of this subject
# thetaInit[sIdx] = sum(yThisSubj)/length(yThisSubj) # proportion
# }
# initsList = list( theta=thetaInit )
# Option 2: Use function that generates random values near MLE:
initsList = function() {
thetaInit = rep(0,Nsubj)
for ( sIdx in 1:Nsubj ) { # for each subject
includeRows = ( s == sIdx ) # identify rows of this subject
yThisSubj = y[includeRows] # extract data of this subject
resampledY = sample( yThisSubj , replace=TRUE ) # resample
thetaInit[sIdx] = sum(resampledY)/length(resampledY)
}
thetaInit = 0.001+0.998*thetaInit # keep away from 0,1
return( list( theta=thetaInit ) )
}
247:卵の名無しさん
17/11/11 16:51:08.92 AQmX3Wf1.net
#-----------------------------------------------------------------------------
# RUN THE CHAINS
parameters = c( "theta") # The parameters to be monitored
adaptSteps = 500 # Number of steps to adapt the samplers
burnInSteps = 500 # Number of steps to burn-in the chains
nChains = 4 # nChains should be 2 or more for diagnostics
thinSteps = 1
nIter = ceiling( ( numSavedSteps * thinSteps ) / nChains )
# Create, initialize, and adapt the model:
jagsModel = jags.model( "TEMPmodel.txt" , data=dataList , inits=initsList ,
n.chains=nChains , n.adapt=adaptSteps )
# Burn-in:
cat( "Burning in the MCMC chain...\n" )
update( jagsModel , n.iter=burnInSteps )
# The saved MCMC chain:
cat( "Sampling final MCMC chain...\n" )
codaSamples = coda.samples( jagsModel , variable.names=parameters ,
n.iter=nIter , thin=thinSteps )
# resulting codaSamples object has these indices:
# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]
if ( !is.null(saveName) ) {
save( codaSamples , file=paste(saveName,"Mcmc.Rdata",sep="") )
}
return( codaSamples )
}
248:卵の名無しさん
17/11/11 19:06:46.40 AQmX3Wf1.net
library("rstan")
library("rstudioapi")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
G1314model='
data{// Golgo13.14.stan
int N13; // 100
int N14; //10
int<lower=0,upper=1> G13[N13];
int<lower=1,upper=2> G14[N14];
}
parameters{
real<lower=0,upper=1> p13;
real<lower=0,upper=1> p14;
}
model {
for(i in 1:N13) G13[i] ~ bernoulli(p13);
for(j in 1:N14) G14[j] ~ bernoulli(p14);
p13 ~ beta(1,1);
p14 ~ beta(1,1);
}
generated quantities{
real d = p13 - p14;
}
'
writeLines( G1314model , con='G1314.stan' )
249:卵の名無しさん
17/11/11 19:15:36.64 AQmX3Wf1.net
N13=100
N14=10
G13=rep(1,N13)
G14=rep(1,N14)
data=list(N13=N13,N14=N14,G13=G13,G14=G14)
G1314.model=stan_model('G1314.stan')
saveRDS(G1314.model,file='G1314.model.rds')
# G1314.model=readRDS('G1314.model.rds')
fitG1314 <- sampling(G1314.model, data=data, seed=1234)
fitG1314
stan_dens(fitG1314,separate_chains = TRUE)
i.imgur.com/aL2PSDM.png
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
p13 0.99 0.00 0.01 0.97 0.99 0.99 1.00 1.00 2500 1
p14 0.92 0.00 0.08 0.72 0.88 0.94 0.97 1.00 2834 1
d 0.07 0.00 0.08 -0.01 0.02 0.05 0.11 0.27 2837 1
lp__ -10.24 0.04 1.16 -13.44 -10.70 -9.88 -9.42 -9.09 979 1
> ms=rstan::extract(fitG1314)
> mean(ms$d<0)
[1] 0.09425
250:卵の名無しさん
17/11/11 19:25:33.55 AQmX3Wf1.net
URLﾘﾝｸ(i.imgur.com)
251:卵の名無しさん
17/11/11 20:30:46.13 AQmX3Wf1.net
Golgo13 vs Golgo14
Let's MCMC with JAGS.
I set the ROPE ( Range of Practical Equivalence ) as -0.01 to 0.01
> smryMCMC( mcmcCoda2 , compVal=NULL ,ropeDiff = c(-0.01,0.01))
Mean Median Mode HDIlow HDIhigh CompVal PcntGtCompVal
Diff 0.07384483 0.05247579 0.01267166 -0.02478474 0.2345431 0 90.18
PcntLtROPE PcntInROPE PcntGtROPE
Diff 3.07 15.66 81.27
It'll be illustrated as follows.
URLﾘﾝｸ(i.imgur.com)
They are about equal to the result with Stan.
> mean(d)
[1] 0.07288324
> median(d)
[1] 0.05019718
> MAP(d)[1]# mode
0.01531325
> HDInterval::hdi(d)
lower upper
-0.0242909 0.2395352
> mean(d>0) # PcntGtCompVal
[1] 0.90575
> mean(d < -0.01) # PcntLtROPE
[1] 0.03475
> mean(-0.01 < d & d < 0.01) PcntInROPE
[1] 0.15975
> mean(d>0.01) # PcntGtROPE
[1] 0.8055
252:卵の名無しさん
17/11/11 23:28:55.02 AQmX3Wf1.net
0/5 vs 1/20
URLﾘﾝｸ(i.imgur.com)
253:卵の名無しさん
17/11/12 16:15:23.53 OKdyWAPj.net
posterior ∝ likelihood * prior
URLﾘﾝｸ(i.imgur.com)
254:卵の名無しさん
17/11/13 10:06:43.47 9LmqAQrY.net
K=12
omega1=0.25
omega2=0.75
a1 = omega1*(K-2)+1
b1 = (1-omega1)*(K-2)+1
a2 = omega2*(K-2)+1
b2 = (1-omega2)*(K-2)+1
curve(dbeta(x,a1,b1))
curve(dbeta(x,a2,b2))
curve(dbeta(x,a1,b1)+dbeta(x,a2,b2))
pdf <- function(theta) (dbeta(theta,a1,b1)+dbeta(theta,a2,b2))/2
n=31
d=0.02
h=function(theta,omega){ # some omega must equal to omega1 or omega2
if(omega>=omega1-d & omega<=omega1+d){
a1 = omega*(K-2)+1
b1 = (1-omega)*(K-2)+1
return(dbeta(theta,a1,b1))
}
if(omega>=omega2-d & omega<=omega2+d){
a2 = omega*(K-2)+1
b2 = (1-omega)*(K-2)+1
return(dbeta(theta,a2,b2))
}
return(0)
}
255:卵の名無しさん
17/11/13 10:07:02.25 9LmqAQrY.net
prior=matrix(NA,n,n) # with outer ERROR!
for(i in 1:n){
for(j in 1:n){
prior[i,j]=h(theta[i],omega[j])
}
}
image(theta,omega,prior,col = heat.colors(12,0.5))
contour(theta,omega,prior,add=TRUE)
persp(theta,omega,prior,col='skyblue',theta = 30)
library(rgl)
persp3d(theta,omega,prior,col='skyblue')
s=6 ; f=3
h2=function(theta,omega) theta^s*(1-theta)^f
likelihood=outer(theta,omega,h2)
image(theta,omega,likelihood,col = terrain.colors(12,0.5))
contour(theta,omega,likelihood,add=TRUE)
persp(theta,omega,likelihood,col='wheat',theta = 30)
library(rgl)
persp3d(theta,omega,likelihood,col='wheat')
posterior=prior*likelihood
image(theta,omega,posterior,col = topo.colors(12,0.5))
contour(theta,omega,posterior,add=TRUE)
persp(theta,omega,posterior,col='maroon',theta = 30)
library(rgl)
persp3d(theta,omega,posterior,col='maroon')
256:卵の名無しさん
17/11/13 12:14:38.18 9LmqAQrY.net
ド底辺シリツ医大は裏口入学と学力で入った例外入学がいるとする。
高卒レベルの基礎学力テストをしたところ裏口入学は不合格率の最頻値が0.75、例外者のそれは0.25であった。
いずれの分布も形状母数和が12のベータ分布に従っていた。
ド底辺シリツ医大でテストしたところ6人が不合格、3人が合格であったとき、ド底辺シリツ医大の裏口入学者の割合を推測せよ。
257:卵の名無しさん
17/11/13 20:42:17.05 JnL2ZVyT.net
omega1=0.25
omega2=0.75
K=12
modelString='
model {
f o r(ii n1 : N){
y[i] ~ dbern( theta )
}
theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )
omega[1] <-
omega1
omega[2] <- omega2
kappa <- K
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1 - p
p ~ dunif(0,1)
}
'
dataList=list(omega1=omega1,omega2=omega2,K=K)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('omega','kappa','theta'), n.iter=10000)
summary(codaSamples)
plot(codaSamples)
258:卵の名無しさん
17/11/13 20:51:13.12 JnL2ZVyT.net
# in debug
y= c(1,1,1,1,1,1,0,0,0)
N=length(y)
omega1=0.25
omega2=0.75
K=12
modelString='
model {
for(i n1:N){
y[i] ~ dbern( theta )
}
theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )
omega[1] <- omega1
omega[2] <- omega2
kappa <- K
m ~ dcat(mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1 - p
p ~ dunif(0,1)
}
'
dataList=list(y=y,N=N, omega1=omega1,omega2=omega2,K=K)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('m,'theta'), n.iter=10000)
summary(codaSamples)
plot(codaSamples)
259:卵の名無しさん
17/11/13 21:18:42.89 JnL2ZVyT.net
# in debug
y= c(1,1,1,1,1,1,0,0,0)
N=length(y)
omega1=0.25
omega2=0.75
kappa1=12
kappa2=12
modelString='
model {
for (i in 1:N){
y[i] ~ dbern( theta )
}
theta <- equals(m,1)*theta1 + equals(m,2)*theta2
theta1 ~ dbeta( omega1*(kappa1-2)+1 , (1-omega1)*(kappa1-2)+1 )
theta2 ~ dbeta( omega2*(kappa2-2)+1 , (1-omega2)*(kappa2-2)+1 )
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1-p
p ~ dunif(0,1)
}
'
dataList=list(y=y,N=N, omega1=omega1,omega2=omega2, kappa1=kappa1, kappa2=kappa2)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('m,'theta'), n.iter=10000)
summary(codaSamples)
plot(codaSamples)
260:卵の名無しさん
17/11/13 21:51:39.83 JnL2ZVyT.net
P(Uraguchi|Wrong) = P(Wrong|Uraguchi)*P(Uraguchi)/(P(Wrong|Uraguchi)*P(Uraguchi)+P(Wrong|Square)*P(Square))
Uraguchi: Buying the way
261:to Do-Teihen, unfair matriculation Wrong: Writing Wrong English Square: fair matriculation P(Square) = 1 - P(Uraguchi) P(Wrong|Uraguchi) = 1 P(Wrong|Square)=0.001 P(Uraguchi) ~ dunif(0,1) P(U|W)=1*P(U)/(1*P(U)+0.001*(1-P(U))) modelString=' puw <- 1 * pu /(1 * pu + 0.001*(1-pu)) pu ~ dunif(0,1) '
262:卵の名無しさん
17/11/14 07:16:26.27 kN15uX//.net
>>249
> js=as.matrix(codaSamples)
> boxplot(theta~m)
> tapply(theta,m,length)
1 2
8778 41222
> sum(m==2)/length(m)
[1] 0.82444
263:卵の名無しさん
17/11/14 11:23:12.95 /SJKjWLk.net
require(rjags)
JmodelString='
model {
puw = 1 * pu /(1 * pu + 0.001*(1-pu))
pu ~ dunif(0,1)
}
'
writeLines( JmodelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt")
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('pu','puw'), n.iter=100000)
summary(codaSamples)
js=as.matrix(codaSamples)
hist(js[,'puw'],xlim=c(0.9,1),freq=FALSE,breaks=1000,col='red',
main='P(Uraguchi|Wrong_English)',xlab='Probability')
HDInterval::hdi(js[,'puw'])
264:卵の名無しさん
17/11/14 13:59:12.73 /SJKjWLk.net
# Stan for Uraguchi
modelString='
parameters{
real<lower=0,upper=1> pu;
}
transformed parameters{
real puw = 1 * pu /(1 * pu + 0.001*(1-pu));
}
model{
pu ~ uniform(0,1);
}
ura.model<-stan_model(model_code = modelString)
fit.ura <- rstan::sampling(ura.model,seed=123,iter=10000,warmup=5000)
print(fit.ura,digits=4)
ms=rstan::extract(fit.ura)
round(HDInterval::hdi(ms$puw),4)
265:卵の名無しさん
17/11/14 18:44:34.87 Njwo3HI0.net
ド底辺シリツ医大は裏口入学と学力で入った例外入学がいるとする。
高卒レベルの基礎学力テストをしたところ裏口入学は不合格率の最頻値が0.75、例外者のそれは0.25であった。
いずれの分布も形状母数和が12のベータ分布に従っていた。
あるド底辺シリツ医大でテストしたところ6人が不合格、3人が合格であったとき、このド底辺シリツ医大の裏口入学者の割合を推測せよ。
Suppose there are Uraguchi and non-Uraguchi(irregular) enrollees in the DoTeihen medical school,
they get the achievement test for high school students.
The failing rate of Uraguchi is known to distribute in β distributuion with its mode value ω = 0.75 and sum of shape parameters κ = 12.
The failing rate of irregulars is in β distribution with ω = 0.25 and κ = 12.
At a DoTeihen medical school, 9 alimni had the test, and 6 failed and 3 passed.
Infer the proportion of Uraguchi in this DoTeihen.
266:卵の名無しさん
17/11/15 10:17:36.22 4qVYF7j+.net
y= c(1,1,1,1,1,1,1,1,1)
N=length(y)
omega1=0.75
omega2=0.25
kappa1=12
kappa2=12
modelString='
model {
for (i in 1:N){
y[i] ~ dbern( theta )
}
theta <- equals(m,1)*theta1 + equals(m,2)*theta2
theta1 ~ dbeta( omega1*(kappa1-2)+1 , (1-omega1)*(kappa1-2)+1 )
theta2 ~ dbeta( omega2*(kappa2-2)+1 , (1-omega2)*(kappa2-2)+1 )
m ~ dcat( mPriorProb[] )
mPriorProb[1] <- p
mPriorProb[2] <- 1-p
p ~ dunif(0,1)
}
'
dataList=list(y=y,N=N, omega1=omega1,omega2=omega2, kappa1=kappa1, kappa2=kappa2)
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('m','theta','theta1','theta2'), n.iter=50000)
summary(codaSamples)
js=as.matrix(codaSamples)
tapply(js[,'theta'],js[,'m'],length)
sum(js[,'m']==1)/nrow(js)
267:卵の名無しさん
17/11/15 14:20:17.72 RUmagzE5.net
In Bayesian data analysis, evidence is the marginal likelihood （Integrate P(D|Θ)P(Θ)dΘ） which MCMC cannot yield.
268:卵の名無しさん
17/11/17 12:56:37.51 nii6SWM6.net
# randam numbers by inverse culmutive density function
RandICDF <- function(ICDF,PDF,...){
U=runif(10000)
rand=ICDF(U,...)
hist(rand,freq=FALSE,breaks=30,
col=sample(colors(),2),main='')
curve(PDF(x,...),add=TRUE,lwd=2)
invisible(rand)
}
par(mfrow=c(2,2))
RandICDF(qnorm,dnorm)
RandICDF(qbeta,dbeta,shape1=3.5,shape2=8.5)
RandICDF(qbeta,dbeta,shape1=0.1,shape2=0.1)
RandICDF(qgamma,dgamma,shape=3,rate=1)
URLﾘﾝｸ(i.imgur.com)
269:卵の名無しさん
17/11/17 13:24:40.98 nii6SWM6.net
# randam numbers following PDF by Jhon von Neuman's method
vonNeumann <- function(PDF,xmin=0,xmax=1){
N=10000
ymax=max(PDF(seq(xmin,xmax,length=N+1)))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,breaks=30,col=sample(colors(),2),main='')
curve(PDF,add=TRUE,lwd=2)
invisible(Rand)
}
par(mfrow=c(2,2))
vonNeumann(function(x)dbeta(x,3.5,8.5))
vonNeumann(function(x)dgamma(x,3,1),0,10)
vonNeumann(dnorm,xmin=-3,xmax=3)
p=runif(1)
f=function(x) p*dnorm(x,-3,1)+(1-p)*dnorm(x,3,3)
vonNeumann(f,xmin=-10,xmax=10)
URLﾘﾝｸ(i.imgur.com)
270:卵の名無しさん
17/11/17 16:06:01.69 nii6SWM6.net
Golgo Script will be reduced as follows:
N shoots with z hits
Nz <- function(N,z,...){
curve(dbeta(x,1+z,1+N-z),xlab='Probability of Hit',
ylab=paste('Probabilty of',z,'
271:Hits out of ',N,'Shoots'),...) hdi=HDInterval::hdi(qbeta,shape1=1+z,shape2=1+N-z) print(hdi,digits=4) }
272:卵の名無しさん
17/11/18 10:52:47.33 V/eIhsOZ.net
modelString =
'model {
z ~ dbin(0.5,N)
N ~ dpois(lambda)
p=z/N
}
'
writeLines(modelString , con="TEMPmodel.txt" )
f <- function (lambda){
dataList=list(lambda=lambda)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('p'), n.iter=10000)
js=as.matrix(codaSamples)
jsp=js[,'p']
mean(jsp<=7/24)
}
xx=7:50
yy=sapply(xx,f)
plot(xx,yy)
273:卵の名無しさん
17/11/18 16:59:45.52 FiDbN6uZ.net
graphics.off()
p=binom.test(7,24,alt='less')$p.value
1-(1-p)^2 # 0.06289339
z1=7
N1=24
N2=24
plot(0:N1/N1,dbinom(0:N1,N1,0.5),type='h',lwd=5,col='skyblue',ann=FALSE,ylim=c(0,0.25))
points(z1/N1,0,pch='+',cex=2)
binom.test(z1,N1,alt='less')$p.value
0:N1/N1 < z1/N1
sum(0:N1/N1 <= z1/N1) # 8個0:7
sum(dbinom(0:z1,N1,0.5))
# A∪B == A + B - A∩B 0.06289339
sum(dbinom(0:7,N1,0.5))+sum(dbinom(0:7,N2,0.5)) -sum(outer(dbinom(0:7,N1,0.5),dbinom(0:7,N2,0.5),'*'))
z1=7
N1=24
N2=12
plot(0:N2/N2,dbinom(0:N2,N2,0.5),type='h',lwd=5,col='wheat',ann=FALSE)
lines(0:N1/N1+0.05,dbinom(0:N1,N1,0.5),type='h',lwd=5,col='skyblue',ann=FALSE,ylim=c(0,0.25))
points(z1/N1,0,pch='+',cex=2)
0:N2/N2 < z1/N1
sum(0:N2/N2 <= z1/N1) # 4個0:3
# A∪B == A + B - A∩B 0.1026226
sum(dbinom(0:7,N1,0.5))+sum(dbinom(0:3,N2,0.5)) - sum(outer(dbinom(0:z1,N1,0.5),dbinom(0:3,N2,0.5)))
274:卵の名無しさん
17/11/20 19:17:52.08 xJug4kDO.net
#
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=FALSE){
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
cdf <- function(x) integrate(pdf,xmin,x)$value/AUC
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]))
if(Print){
pp=seq(0,1,length=nxx)
plot(pp,sapply(pp,ICDF),type='l',lwd=2,xlab='p',ylab='x')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
}
invisible(ICDF)
}
pdf2hdi(function(x)dbeta(x,0.001,0.001)*x^7*(1-x)^17)
pdf2hdi(dnorm,-10,10,cre=0.95)
275:卵の名無しさん
17/11/20 20:38:04.94 xJug4kDO.net
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=4)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='l',lwd=2,xlab='x',ylab='Density')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='l',lwd=2,xlab='x',ylab='Probability')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',lwd=2,xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}
276:卵の名無しさん
17/11/20 21:23:09.45 xJug4kDO.net
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)
}
277:卵の名無しさん
17/11/20 22:09:25.75 xJug4kDO.net
According to this posting, ｽﾚﾘﾝｸ(doctor板:529番) as many as 15 freshmen flunk.
Let's assume there are 120 students in one class and the number of flunker is distributed as poisson distribution.
In what range will students are expected to flunk next year? Calculate the number with 99% confidence interval.
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE,FUN=FALSE){
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))
}
if(FUN) invisible(ICDF)
invisible(hdi)
}
278:卵の名無しさん
17/11/22 10:27:31.38 WHgLvp/3.net
ab2mv<-function(a,b){
m<-a/(a+b)
v<-m*(1-m)/(a+b+1)
mv<-c(m,v)
return(mv)
}
mv2ab<-function(m,v){
a=(-m^3+m^2-m*v)/v
b=(m^3-2*m^2+m*v+m-v)/v
ab<-c(a,b)
return(ab)
}
HDCI <- function(PMF,cl=0.95){ # Highest Density Confidence Interval
PDF=PMF/sum(PMF)
rsPDF=rev(sort(PDF))
min.density=rsPDF[min(which(cumsum(rsPDF)>=cl))]
index=which(PDF>=min.density)
data.frame(lower.idx=round(min(index)),upper.idx=round(max(index)),actual.CI=sum(PDF[index]))
}
279:卵の名無しさん
17/11/22 12:17:28.03 WHgLvp/3.net
Scale <- function(x) (x-mean(x))/sqrt(sum((x-mean(x))^2)/(length(x)-1))
m=50;s=10
rn
280:=rnorm(1000,m,s) mean(rn) sd(rn) y=Scale(rn)*s+m mean(y) sd(y)
281:卵の名無しさん
17/11/22 19:57:10.44 WHgLvp/3.net
URLﾘﾝｸ(youtu.be)
282:卵の名無しさん
17/11/23 15:22:14.60 Ue5tZuwc.net
par(mfrow=c(2,2))
theta=0.5
NN=1000
N=1:NN
flip=numeric()
for(i in N){
flip=append(flip,rbinom(1,1,theta))
}
z=cumsum(flip)
z_N=z/N
plot(z_N,type='l',ylim=c(0,1))
pv=numeric()
for(i in N){
pv[i]=binom.test(z[i],N[i],theta)$p.value
}
plot(pv,type='l')
abline(h=0.05,col='blue',ylim=c(0,1))
bf=numeric()
for(i in N){
bf[i]=beta(z[i]+1,N[i]-z[i]+1)/beta(1,1)/(theta^z[i]*(1-theta)^(N[i]-z[i]))
}
plot(log(bf),type='l')
abline(h=log(3),lty=3)
abline(h=log(1/3),lty=3)
283:卵の名無しさん
17/11/23 15:22:43.97 Ue5tZuwc.net
hdi=NULL
for(i in N){
y=flip[1:i]
s=rep(1,i)
data=data.frame(y,s)
Ntotal=i
Nsubj=1
dataList=list(y=y,s=s,Ntotal=Ntotal,Nsubj=Nsubj)
js=genMCMC2(data)
hdi=rbind(hdi,HDIofMCMC(as.matrix(js)))
}
saveRDS(hdi,'hdi_sequential')
hdi=readRDS('hdi_sequential')
plot(hdi[,1],type='l',ylim=c(0,1),main='95% HDI')
lines(hdi[,2])
abline(h=0.5,col=4)
plot(apply(hdi,1,diff),type='l',main='HDI width')
284:卵の名無しさん
17/11/25 12:10:33.25 3kOACkIe.net
Metropolis Algo
source('DBDA2E-utilities.R')
Bern <- function(x) rbinom(1,1,x)
Metro <- function(
SD=0.02,
.z=14,
.N=20,
a=1,
b=1,
k=50000,
Print=TRUE){
theta=numeric()
theta[1]=0.01
likely <- function(x,z=.z,N=.N) x^z*(1-x)^(N-z)
for(i in 1:(k-1)){
delta=rnorm(1,0,SD)
theta_p=theta[i]+delta
p=min(1,likely(theta_p)*dbeta(theta_p,a,b)/
(likely(theta[i])*dbeta(theta[i],a,b)))
theta[i+1]=theta[i]+Bern(p)*delta
}
if(Print){
# plot(theta,1:k,type='l')
plotPost(theta,cenTend = 'mean',cex.lab = 1)
plot(theta[1:100],1:100,type='l',xlim=c(0,1))
plot(theta[(k-100):k],(k-100):k,type='l',xlim=c(0,1))}
invisible(theta)
}
285:卵の名無しさん
17/11/26 19:08:27.54 11LwtYgG.net
# JAGS for proportion
graphics.off()
rm(list=ls())
zi=100; Ni=100; zj=10; Nj=10
(y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj)))
(s=as.numeric(factor(c(rep('D',Ni),rep('U',Nj)))))
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
a=1 ; b=1 # prior : beta(a,b)
# JAGS model
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( sIdx in 1:Nsubj ) {
theta[sIdx] ~ dbeta(", a,',' , b," )
}
}
")
# close modelString
286:卵の名無しさん
17/11/26 19:09:02.05 11LwtYgG.net
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
#jags.model(file, data, inits,n.chains = 1, n.adapt=1000, quiet=FALSE)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
#coda.samples(model, variable.names, n.iter, thin = 1, na.rm=TRUE, ...)
summary(codaSamples)
plot(codaSamples,col=sample(colors()))
mcmcMat=as.matrix(codaSamples)
par(mfrow=c(2,2))
hist(mcmcMat[,1],freq=FALSE,xlim=c(0,1),
col=sample(colours(),1))
hist(mcmcMat[,1]-mcmcMat[,2],freq=FALSE,xlim=c(-1,1),
col=sample(colours(),1))
plot(mcmcMat[,1],mcmcMat[,2], col=rgb(0.01,0.01,0.3,0.25))
hist(mcmcMat[,2],freq=FALSE, xlim=c(0,1),
col=sample(colours(),1))
287:卵の名無しさん
17/11/26 19:09:08.75 11LwtYgG.net
source("Kruschke_tools.R") # for genMCMC2
(myData=data.frame(y,s))
mcmcCoda = genMCMC2( data=myData , numSavedSteps=10000,a=1,b=1)
# mcmcCoda = genMCMC( data=myData , numSavedSteps=10000)
# genMCMC
# Display diagnostics of chain, for specified parameter:
source('Jags-Ydich-XnomSsubj-MbernBeta.R')
diagMCMC( mcmcCoda , parName="theta[1]" )
diagMCMC( mcmcCoda , parName="theta[2]" )
# Display numerical summary statistics of chain:
smryMCMC( mcmcCoda , compVal=NULL,ropeDiff = c(-0.025,0.025))
#function( codaSamples , compVal=NULL , rope=NULL , saveName=NULL )
summary(mcmcCoda)
# Display graphical posterior information:
plotMCMC( mcmcCoda , data=myData , compVal=NULL,ropeDiff = c(-0.025,0.025))
# function( codaSamples , data , compVal=NULL , rope=NULL ,
# saveName=NULL , showCurve=FALSE , saveType="jpg" )
plot(mcmcCoda)
plotMCMC2( mcmcCoda , data=myData , compVal=NULL, showCurve=FALSE,
.credMass = 0.95,ropeDiff = c(-0.025,0.025),cenTend='mean')
(summry=smryMCMC( mcmcCoda, compVal=NULL, ropeDiff = c(-0.025,0.025)))
print(summry[3,],digits=2)
prop.test(c(zi,zj),c(Ni,Nj))$p.value
fisher.test(matrix(c(zi,Ni-zi,zj,Nj-zj),2))$p.value
288:卵の名無しさん
17/11/28 18:49:23.50 QFCizNPN.net
logistic <- function(x,gain,threshol
289:d){ 1/(1 + exp(-gain*(x-threshold))) } b2gt <- function(b0,b1,b2){ gain=sqrt(b1^2+b2^2) threshold=-b0/gain return(gain=gain,threshold=threshold) }
290:卵の名無しさん
17/11/28 18:52:07.55 QFCizNPN.net
I could not locate a good site to explain normalizationn for logististic regression,
but with the examples depicted in the textbook I have finally got understood.
This is it.
URLﾘﾝｸ(i.imgur.com)
core portion of its code :
ｽﾚﾘﾝｸ(hosp板:275番)
291:卵の名無しさん
17/11/28 20:13:41.22 QFCizNPN.net
# z=b0+巴k*xk
# p = logistic(z) = 1/(1+e^-z)
# 1-p = (1+e^-z)/(1+e^-z) - 1/(1+e^-z) = e^-z/(1+e^-z)
# p/(1-p) = 1 /e^-z = e^z
# logit(p) = log(p/(1-p))= log(e^z) = z
# logit(logistic(z)) = z
292:卵の名無しさん
17/11/30 16:34:37.18 TyAFrmPC.net
dataList=list(y=y,Ntotal=length(y),meanY=mean(y),sdY=sd(y))
modelString = '
model {
for ( i in 1:Ntotal ) {
y[i] ~ dt( mu , 1/sigma^2 , nu )
}
mu ~ dnorm( meanY , 1/(100*sdY)^2 )
sigma ~ dunif( sdY/1000 , sdY*1000 )
nu ~ dexp(1/30.0)
}
'
writeLines(modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable=c('sigma','mu','nu'), n.iter=10000)
plot(codaSamples,col=sample(colours(),1))
js=as.matrix(codaSamples)
head(js)
293:卵の名無しさん
17/11/30 16:35:08.68 TyAFrmPC.net
# Y = aX + b , X ~ dt, a:scale parameter, b:location parameter
dt_ls <- function(x, df, mu, a) 1/a * dt((x - mu)/a, df)
pt_ls <- function(x, df, mu, a) pt((x - mu)/a, df)
qt_ls <- function(prob, df, mu, a) qt(prob, df)*a + mu
rt_ls <- function(n, df, mu, a) rt(n,df)*a + mu
par(mfrow=c(1,1))
hist(y,breaks=20,col='skyblue',freq=FALSE,xlim=c(30,220),main='')
N=63 #length(y)
for(i in sample(1:nrow(js),N)){
curve(dt_ls(x,js[i,'nu'],js[i,'mu'],js[i,'sigma']),add=TRUE,
lty=1,col=rgb(.01,.01,.01,.1))
}
294:卵の名無しさん
17/12/01 13:08:02.18 UfpWtEOZ.net
dT <- function(x, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
dt((x - mu)/s, nu)/s
}
pT <- function(x, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
pt((x - mu)/a, nu)
}
qT <- function(prob, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
qt(prob, nu)*s + mu
}
rT <- function(n, nu, mu, sd){
s=sd*sqrt((nu-2)/nu)
rt(n,nu)*s + mu
}
295:卵の名無しさん
17/12/02 06:06:16.06 SDqtqHE2.net
男性 28.2％
女性 9.0％
男女計 18.2％
URLﾘﾝｸ(www.jti.co.jp)
P(s)=.182
P(s|f)=.090
P(s|m)=.282
P(s)=P(s|f)P(f)+P(s|m)P(m)
P(m)=1-P(f)
から
P(f)=(P(s)-P(s|m))/(P(s|f)-P(s|m))
ベイズの公式
P(f|s)=P(s|f)P(f)/(P(s|m)P(m)+P(s|f)P(f))
P(s|m)P(m)+P(s|f)P(f)=P(s)
.090*(.182-.282)/(.090-.282)/0.182=0.2575549
296:卵の名無しさん
17/12/02 11:15:07.32 ZaK9sW49.net
# 問.
# 患者が煙草を忘れて行ったとする。
# 忘れて行った人物が女性である確率を以下のデータから計算せよ。
#
# 喫煙率
# 男性 28.2％
# 女性 9.0％
# 男女計 18.2％
P(s|m) = 0.282
P(s|f) = 0.090
P(s) = P(s|m)P(m)+P(s|f)P(f)
= P(s|m)(1-P(f)) + P(s|f)P(f)
= 0.182
P(f) = (P(s) - P(s|m))/(P(s|f) - P(s|m))
= (0.182 - 0.282)/(0.090 - 0.282)
= 0.5208333
P(f|s) = P(s|f)P(f)/P(s)
= 0.090*0.5208333/0.182
= 0.2575549
#
LR = P(s|f)|P(s|m)=0.090/0.282=0.3191489
prior.odds(f)=P(f)/(1-P(f))=0.5208333/(1-0.5208333)=1.086956
post.odds(f|s)= prior.odds(f)*LH=1.086956*0.3191489=0.3469008
P(f|s)=post.odds(f|s)/(1+post.odds(f|s))=0.3469008/(1+0.3469008)
= 0.2575548
297:卵の名無しさん
17/12/03 06:11:30.00 B6LMarvh.net
１次方程式もできないド底辺特殊シリツ医大卒の記録
URLﾘﾝｸ(imagizer.imageshack.com)
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ？
裏口入学以外に説明がつく？
中学生でも解ける一次方程式の問題だろ。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか？
だとさ。
URLﾘﾝｸ(imagizer.imageshack.com)
0.2575549
と答を書いてやったら
>単位も書かずに答えだとか…
ド底辺シリツ医大では確率に単位があるらしいぞｗｗｗ
何でこんな馬鹿が大学に入れるわけ？
裏口入学以外に説明がつく？
URLﾘﾝｸ(imagizer.imageshack.com)
298:卵の名無しさん
17/12/03 08:22:35.16 Egt6Q5KK.net
１次方程式もできないド底辺特殊シリツ医大卒の記録
URLﾘﾝｸ(imagizer.imageshack.com)
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ？
裏口入学以外に説明がつく？
中学生でも解ける一次方程式の問題だろ。
シリツ医大には二次方程式が解けないやつがいると言ってた えなりかずき もビックリだろね。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか？
だとさ。
URLﾘﾝｸ(imagizer.imageshack.com)
求める確率を
0.2575549
と答を書いてやったら
>単位も書かずに答えだとか…
ド底辺シリツ医大では確率に単位があるらしいぞｗｗｗ
何でこんな馬鹿が大学に入れるわけ？
裏口入学以外に説明がつく？
URLﾘﾝｸ(imagizer.imageshack.com)
299:卵の名無しさん
17/12/03 10:56:29.84 qW8l0b6t.net
# ある仮想の難治疾患患者25人従来薬を投与して3人治癒した。
# 新薬が登場して3人に投与したところ治癒した人はいなかった。
# この新薬を継続して使う価値があるかどうか検討せよ。
別バージョン
# 巨乳女子大で25人に声をかけたら3人が誘いにのった。
# 桃尻女子大で3人に声をかけたら誰も誘いにのらなかった。
# どちらが口説きやすいか検討せよ。
JAGSでMCMCして治癒率の確率密度関数を描くとこうなる。
URLﾘﾝｸ(i.imgur.com)
治癒率差の不偏推定量は
> mean(dif)
[1] -0.05136971
54％が負
> c(mean(dif<0),mean(dif>0))
[1] 0.5395 0.4605
５％幅の違いは同等扱いにすると
> c(mean(dif<ROPE[1]),mean(ROPE[1]<dif & dif<ROPE[2]), mean(dif>ROPE[2]))
[1] 0.4834 0.1236 0.3930
と計算できる。
98％HDIは
> HDInterval::hdi(dif)
lower upper
-0.4247349 0.2535357
と０を挟む。
RのパッケージBESTを改造して、治癒率の差の関数密度をかくと
URLﾘﾝｸ(i.imgur.com)
ゆえに新薬は無効とはいえないだけでなく、不偏推定量から従来薬を凌駕する可能性が54％ある。