17/09/27 16:26:42.51 By0OihfH.net
N=100 ; n=5 ; r=0 ; m=10
atari.Q <- function(N,n,r,k=10^5){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
print(c(mean = sum(xx*pdf),
mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1))
cat('\n')
print(quantile(SS,c(0.025,0.5,0.975)))
invisible(SS)
}
lottery.mci <- function(N,m,k=10^5){
pmf <- function(S,.N=N,.m=m){
choose(.N-.m,S-1)/choose(.N,S)
}
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:N,N*k,replace=TRUE,prob=pdf)
print(c(mean=E,mod=which.max(pdf)))
cat('\n')
print(quantile(SSS,c(.025,.975)))
invisible(SSS)
}
atari.Q(N,n,r)
lottery.mci(N,m)
101:卵の名無しさん
17/09/27 16:27:31.48 By0OihfH.net
>>97
シリツ医大進学予備校"どていへん予備校"では５年連続不合格、
別のシリツ医大進学予備校"うらぐち予備校"では１０年めで初めて合格者がでたという実績があるとき、
どていへん予備校 と うらぐち予備校ではどちらが合格可能性が高いと言えるか？
102:innuendo
17/09/27 16:37:18.00 By0OihfH.net
巨乳女子大100人と桃尻女子大100人のどちらが誘いに応じやすいか検討してみる。
巨乳女子大では10人に声をかけて2人が誘いに応じた。
桃尻女子大では無作為順に声をかけていったら5人めが誘いに応じた。
どちらが誘いに応じやすいといえるか？
URLﾘﾝｸ(i.imgur.com)
> lottery.mci(100,5)
mean mod
28.14286 21.00000
2.5% 97.5%
4 63
> atari.Q(100,10,2)
mean mode
24.5 20.0
2.5% 50% 97.5%
6 23 50
>
103:卵の名無しさん
17/09/27 17:32:04.45 By0OihfH.net
N=100 ; n=5 ; r=0 ; m=10
atari.Q <- function(N,n,r,k=10^5){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
print(c(mean = sum(xx*pdf),
mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1))
cat('\n')
print(quantile(SS,c(0.025,0.5,0.975)))
invisible(SS)
}
lottery.mci <- function(N,m,k=10^5){
pmf <- function(S,.N=N,.m=m){
choose(.N-.m,S-1)/choose(.N,S)
}
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:N,N*k,replace=TRUE,prob=pdf)
print(c(mean=E,mod=which.max(pdf)-1))
cat('\n')
print(quantile(SSS,c(.025,.975)))
invisible(SSS)
104:卵の名無しさん
17/09/27 17:40:20.17 By0OihfH.net
N=100 ; n=10 ; r=2 ; m=5 # バグ修正版
atari.Q <- function(N,n,r,k=10^5){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
# print(quantile(SS,c(0.025,0.5,0.975)))
# cat('\n')
# print(c(mean = sum(xx*pdf),
# mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1))
invisible(SS)
}
lottery.mci <- function(N,m,k=10^5){
pmf <- function(S,.N=N,.r=m){
choose(.N-.r,S-1)/choose(.N,S)
}
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:N,N*k,replace=TRUE,prob=pdf)
UL=quantile(SSS,c(.025,.975))
# print(c(E,UL))
invisible(SSS)
}
Brown=sum(choose(0:N,r)*choose(N-0:N,n-r)/choose(N,n))
curve(choose(x,r)*choose(N-x,n-r)/choose(N,n)/Brown,0,100,col='brown',xlab='尻軽数',ylab='')
Pink=sum(choose(N-m,0:N-1)/choose(N,0:N))
curve(choose(N-m,x-1)/choose(N,x)/Pink,0,100, add=TRUE,col='pink',lwd=2)
legend('topright',bty='n',legend=c('巨乳女子','桃尻女子'),lwd=1:2,col=c('brown','pink'))
105:innuendo
17/09/27 17:59:40.63 By0OihfH.net
URLﾘﾝｸ(i.imgur.com)
平均値と信頼区間のどちらを選択するかだな。
106:卵の名無しさん
17/09/27 20:59:53.74 By0OihfH.net
options(scipen = 100)
lottery.s <- function(m,N=100){
pmf <- function(S,.N=N,.m=m){
choose(.N-.m,S-1)/choose(.N,S)
}
bb=0:N # numbers of possible bitches
pdf=pmf(bb)/sum(pmf(bb))
SS=sample(0:N,replace=TRUE,prob=pdf) # sample as the same length
return(SS)
}
lottery.s(8)
foo <- function(y,al='l'){ # one-sided mo > ky
f <- function(x){
mo=lottery.s(5)
ky=lottery.s(x)
t.test(ky,mo,alt=al)$p.value
}
welch=replicate(10^2,f(y))
mean(welch)
}
yy=5:10
pp=sapply(yy,foo)
plot(yy,pp)
abline(h=0.05,lty=3)
107:卵の名無しさん
17/09/27 21:27:00.63 By0OihfH.net
foo <- function(y,fn=t.test,al='l'){ # one-sided mo > ky
f <- function(x){
mo=lottery.s(5)
ky=lottery.s(x)
fn(ky,mo,alt=al)$p.value
}
p.values=replicate(10^2,f(y))
mean(p.values)
}
yy=5:10
# どの検定でも結果は変わらず
pp=sapply(yy,function(x) foo(x,fn=t.test))
pp=sapply(yy,function(x) foo(x,fn=lawstat::brunner.munzel.test))
pp=sapply(yy,function(x) foo(x,fn=wilcox.test))
PP=sapply(yy,function(x) foo(x,fn=perm::permTS))
plot(yy,pp)
abline(h=0.05,lty=3)
108:卵の名無しさん
17/09/27 23:32:37.60 By0OihfH.net
> lottery.e <- function(N,m){
+ pmf <- function(S,.N=N,.m=m){
+ choose(.N-.m,S-1)/choose(.N,S)
+ }
+ ss=0:N
+ pdf=pmf(ss)/sum(pmf(ss))
+ sum(ss*pdf)
+ }
>
> lottery.e(100,10)
[1] 16
> lottery.e(50,5)
[1] 13.85714
>
> A=100;a=10
> B=50;b=5
>
> prop.test(c(lottery.e(A,a),lottery.e(B,b)),c(A,B))
2-sample test for equality of proportions with continuity correction
data: c(lottery.e(A, a), lottery.e(B, b)) out of c(A, B)
X-squared = 2.1814, df = 1, p-value = 0.1397
alternative hypothesis: two.sided
95 percent confidence interval:
-0.27551116 0.04122545
sample estimates:
prop 1 prop 2
0.1600000 0.2771429
109:卵の名無しさん
17/09/28 00:10:25.54 bltMOtIo.net
sillygal <- function(a,b,A,B){
lottery.e <- function(N,m){
pmf <- function(S,.N=N,.m=m) choose(.N-.m,S-1)/choose(.N,S)
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
sum(ss*pdf)
}
prop.test(c(lottery.e(A,a),lottery.e(B,b)),c(A,B))
}
sillygal(5,10,50,100)
sillygal (5,7,100,100)
110:卵の名無しさん
17/09/29 12:14:08.07 aSJfCsJ7.net
1 st 2 nd 3 rd 4 th と表示させるための1行スクリプト
th=ifelse(0<r&&r<4,c('st','nd','rd')[which(r==1:3)],'th')
111:卵の名無しさん
17/09/29 19:11:58.46 JC0gA3LG.net
N=100
plot(0:N/N,dbinom(0:N,N,p=(0:N)/N),pch='+',
xlab='治癒確率',ylab='治癒確率通り治癒する確率')
yy=dbinom(0:N,N,p=(0:N)/N)
summary(yy[2:100])
hist(yy[2:100],col='tomato')
112:卵の名無しさん
17/09/30 18:41:04.08 ZzasON2B.net
# 元利均等返済 毎月の返済額は一定
.D0=10^8 # 借入額
.ra=0.01 # 年利
.r=.ra/12 # 月利
.n=120 # 返済月数
L <- function(x,D0=.D0,r=.r,n=.n){ # x:毎月の返済額
D=numeric(n)
D[1]=D0*(1+r) - x
for(i in 1:(n-1)) D[i+1] <- D[i]*(1+r) -x
return(D[n]) # n回返済後の借入金残高を返す
}
# L(.n)=0になる毎月の返済額 x を求める。
(m=uniroot(L,c(.D0/.n,.D0*(1+.n*.r)))$root)
m*.n # 総返済額
##### 元利均等返済 毎月の返済額と総返済額
loan <- function(.D0,.ra,.n){
A0=.D0
r=.ra/12
N=.n
m=A0*r*(1+r)^N/( (1+r)^N- 1)
c(Monthly=floor(m),Total=floor(m)*N)
}
loan(23700000,0.03,60)
loan(10^8,0.01,120)
113:卵の名無しさん
17/09/30 18:41:31.88 ZzasON2B.net
## 元利均等返済 nケ月めの残高
Debt <- function(n,.D0,.ra,.n){
A0=.D0
r=.ra/12
N=.n
s=1+r
A0*(1- (1-s^n)/(1-s^N))
}
Debt(1:60,23700000,0.03,60)
##
Loan <- function(.D0,.ra,.n){
L <- function(x,D0=.D0,r=.ra/12,n=.n){
D=numeric(n)
D[1]=D0*(1+r) - x
for(i in 1:(n-1)) D[i+1] <- D[i]*(1+r) -x
return(D[n])
}
m=uniroot(L,c(.D0/.n,.D0*(1+.n*.ra/12)))$root
c(Monthly=floor(m),Total=floor(m)*.n)
}
114:卵の名無しさん
17/10/01 12:38:58.98 NynEuCM1.net
foo <- function(N){
#N=7 # 元数
I=1:N
.Y=sample(I,N) ; .Y # 1,2,..,Nを.Y[1],Y.[2],..,.Y[N]に置換
tikan <- function(X,Y=.Y){ # 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) # 置換後の配列を返す
}
## 循環したらやめる
z=matrix(1:N,nrow=1)
z[1,]=tikan(I)
for(i in 1:N^2){
if(!prod(z[i,]==I)){
z=rbind(z,tikan(z[i,]))
}
else break
}
print(z)
print(rbind(I,z[1,]))
return(nrow(z))
}
115:卵の名無しさん
17/10/01 12:50:45.81 NynEuCM1.net
junkan <- function(.Y){ # 何回目で巡回したか返す、巡回しなければ０
N=length(.Y) # 元数
I=1:N
#.Y=sample(I,N) ; .Y # 1,2,..,Nを.Y[1],Y.[2],..,.Y[N]に置換
tikan <- function(X,Y=.Y){ # 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) # 置換後の配列を返す
}
## 循環したらやめる
z=matrix(1:N,nrow=1)
z[1,]=tikan(I)
for(i in 1:N^2){
if(!prod(z[i,]==I)){
z=rbind(z,tikan(z[i,]))
}
else break
}
# print(rbind(I,.Y))
# print(z)
return(ifelse(i!=N^2,i,ifelse(prod(z[i,]==I),i,0)))
}
116:卵の名無しさん
17/10/01 12:51:32.13 NynEuCM1.net
#
library(gtools)
N=7
PM=permutations(N,N,v=1:N)
jn=numeric(N)
n=nrow(PM)
for(j in 1:n) jn[j] <- junkan(PM[j,])
table(jn)
which(jn==0)
> table(jn)
jn
1 2 3 4 5 6 7 10 12
1 231 350 840 504 1470 720 504 420
> which(jn==0)
integer(0)
117:卵の名無しさん
17/10/01 17:13:29.79 NynEuCM1.net
test,mt1,st1,mt2,st2,mc1,sc1,mc2,sc2
開眼片足立ち左秒,16.02,7.73,18.99,8.12,18.92,9.84,18.83,9.92
開眼片足立ち右秒,18.48,8.29,21.23,8.14,19.88,9.36,19.75,9.35
左握力kg,18.72,4.22,19.67,4.03,19.83,4.63,19.73,4.68
右握力kg,20.89,4.13,21.7,3.86,21.21,4.37,20.92,4.28
FRcm,22.36,6.08,24.83,5.89,23.51,4.97,22.45,4.67
10m歩行速度秒,6.61,0.75,5.95,0.67,6.53,1.18,6.65,1.32
立位体前屈cm,6.90,10.64,12.33,6.41,9.25,4.90,9.00,5.48
左片足立ち振り回,8.40,4.76,9.93,4.24,8.96,6.28,8.36,5.85
右片足立ち振り回,9.30,5.20,11.16,5.34,9.93,7.38,9.80,7.27
nt=30
nc=30
(d=read.csv('clipboard'))
mt2=d[,'mt2']
st2=d[,'st2']
mc2=d[,'mc2']
sc2=d[,'sc2']
# t検定(生データなし,等分散不問)
Welch.test=function(n1,n2,m1,m2,sd1,sd2){
T=(m1-m2)/sqrt(sd1^2/n1+sd2^2/n2)
df=(sd1^2/n1+sd2^2/n2)^2 / (sd1^4/n1^2/(n1-1)+sd2^4/n2^2/(n2-1))
p.value=2*pt(abs(T),df,lower.tail = FALSE)
return(p.value)
}
cbind(d['test'],Welch.test(nt,nc,mt2,mc2,st2,sc2))
118:卵の名無しさん
17/10/03 10:21:03.41 D0GaF7Tg.net
# 百発百中
Dbinom <- function(s,n,p){
choose(n,s)*p^s*(1-p)^(n-s)
}
Dbinom(100,100,0.95)
Dbinom(100,100,0.99)
s=99;n=100 # 百発99中
curve(Dbinom(s,n,x),ylab='',xlab='p')
AUC=integrate(function(p) Dbinom(s,n,p),0,1)$value
integrate(function(x) x*Dbinom(s,n,x)/AUC,0,1)$value
s=100;n=100 # 百発百中
curve(Dbinom(s,n,x),ylab='',xlab='p',0.9,1)
AUC=integrate(function(p) Dbinom(s,n,p),0,1)$value
integrate(function(x) x*Dbinom(s,n,x)/AUC,0,1)$value # Expected value
curve(Dbinom(s,n,x)/AUC)
curve(Dbinom(s,n,x)/AUC,0.9,1)
integrate(function(x) Dbinom(s,n,x)/AUC, 0,1)
pdf <- function(x) Dbinom(s,n,x)/AUC
p99 <- uniroot(function(x,u0=0.99) integrate(pdf,x,1)$value - u0,c(0,1))$root ; p99
N=10^7
hits=rbinom(N,n,p99)
hist(hits)
mean(hits==100)
119:卵の名無しさん
17/10/03 14:56:19.98 D0GaF7Tg.net
# ｎ発ｎ中 CL=0.99
# PMF(x)=nCn*x^n*(1-x)^(n-n)=x^n
# AUC=integrate(PMF,0,1)=[x^(n+1)/(n+1)]1,0 = 1/(n+1)
# PDF=PMF/AUC=(n+1)x^n
# CDF(u)=integrate(pdf(x),0,u)=[x^(n+1)]0,u = u^(n+1)
# (1-CL)^(1/(n+1)) # lower.limit
# exp(log(1-CL)/(n+1))
#
# Ex=integrate(x*PDF,0,1) = integrate((n+1)*x^(n+1),0,1)
# = (n+1)/(n+2)
#
sniper <- function(n,CL=0.99){
c(Lower=(1-CL)^(1/(n+1)), Expected=(n+1)/(n+2), Upper=1)
}
sniper(100,0.99)
nn=1:100
plot(nn,sapply(nn,function(x) sniper(x,CL=0.99)[1]),
xlab='ｎ発ｎ中',ylab='99%信頼区間下限値',pch=19)
120:卵の名無しさん
17/10/03 14:56:47.23 D0GaF7Tg.net
##
N=10^3
n=100
pn=sniper(n,0.99)[1] ; pn
hits=rbinom(N,n,pn)
hist(hits,freq=FALSE, main='')
mean(hits==n)
k=10^5
M=numeric(k)
foo <-function(){
hits=rbinom(N,n,pn)
mean(hits==n)
}
for(i in 1:k) M[i] <- foo()
summary(M)
hist(M,freq=FALSE,main='',xlab='全発命中割合')
121:卵の名無しさん
17/10/04 11:34:22.68 Fg8f7jBH.net
2-way Contingency Table Analysis URLﾘﾝｸ(statpages.info)
Epi::twoby2
122:卵の名無しさん
17/10/04 13:13:51.18 Fg8f7jBH.net
# logRR±1.96*sqrt(1/a-1/(a+b)+1/c-1/(c+d))
# logOR±1.96*sqrt(1/a+1/b+1/c+1/d)
HRCI <- function(a,b,c,d,CL=0.95){
Z=qnorm(1-(1-CL)/2)
HR=(a/(a+b)) /(c/(c+d))
se=sqrt(1/a-1/(a+b)+1/c-1/(c+d))
lwr=HR*exp(-Z*se)
upr=HR*exp( Z*se)
c(lwr,HR,upr)
}
123:innuendo
17/10/04 20:50:12.19 Fg8f7jBH.net
To calculate M1, the relative risks in each of the six studies were combined using a random effects meta-analysis to give a point estimate of 0.361 for the relative risk with a confidence interval of (0.248, 0.527).
The 95% confidence interval upper bound of 0.527 represents a 47% risk reduction, which translates into a risk increase of about 90% from not being on warfarin (1/0.527 = 1.898)
(i.e., what would be seen if the test drug had no effect). Thus, M1 (in terms of the hazard ratio favoring the control to be ruled out) is 1.898.
124:innuendo
17/10/04 21:03:03.78 Fg8f7jBH.net
The clinical margin M2 representing the largest acceptable level of inferiority was therefore set at 50% of M1.
M2 was calculated on the log hazard scale as 1.378, so that NI would be demonstrated provided the upper bound for the 95% confidence interval for C vs. T < 1.378.
exp(0.5*log(1/0.527))=1.378
125:innuendo
17/10/05 11:09:02.87 56lc+jur.net
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
126:卵の名無しさん
17/10/05 13:06:15.38 8CwGr/wg.net
The historical trials assessed Ch ? P, the effect of the active control compared to placebo.
The NI trial assesses T ? Cn , the effect of the test drug compared to active control. If the outcome for the active control is constant across studies,
then Cn = Ch, and the sum T ? Cn + Ch ? P = T ? P represents the effect of the test drug compared to placebo.
127:卵の名無しさん
17/10/05 15:45:07.34 8CwGr/wg.net
# Non-Inferiority Trial
# Fixed Margin Approach
M1M2 <- function(ub){ # ub :upper border of c.i. by meta-analysis
c(M1=1/ub,M2=exp(-log(ub)/2))
}
M1M2(.527)
## Synthesis Method risk ratio
NIs <- function(rr21,se21,rr10,se10){ # 0:placebo, 1:active control 2:test drug
se=sqrt(se21^2+(se10/2)^2) # log_se=sqrt(log(1/a-1/(a+b)+1/c-1/c+d)))
z=(log(rr21)+log(rr10)/2)/se
p=pnorm(-abs(z))
c(p.value=p,z=z)
}
NIs(1.39,0.216,0.361,0.154)
## fixed margin test
NIf <-function(rr21,se21,rr10,se10){
se=se21+se10/2
z=(log(rr21)+log(rr10)/2)/se
p=pnorm(-abs(z))
c(p.value=p,z=z)
}
NIf(1.39,0.216,0.361,0.154)
128:卵の名無しさん
17/10/07 05:54:48.32 427btajj.net
G15 <- function(.N, .n, r, k=10^3,...){
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=0:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1)
}
summary(SS)
}
G(10000,100,100)
129:卵の名無しさん
17/10/07 06:05:26.39 427btajj.net
G13 <- function(N,n,r){
pp=0:N
f <- function(x) choose(x,r)*choose(N-x,n-r)/choose(N,n)
sum(pp*f(pp)/sum(f(pp)))
}
G13(10000,100,100)
G13(10000,10,10)
G13(10000,1,1)
130:卵の名無しさん
17/10/07 16:41:37.63 jwh+lzlc.net
Go13 <- function(.N, .n, r, k=10^3){
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)
}
summary(SS)
}
Go13(100,1,1)
Go13(1000,10,10)
Go13(10000,100,100)
131:卵の名無しさん
17/10/07 19:50:28.87 427btajj.net
Na=1000 ; Nb=1000
pa=0.1 ; pb=0.01
Wa=100 ; Wb=1000
A=c(rep(1,Na*pa),rep(0,Na-Na*pa)) ; A=sample(A)
B=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb)) ; B=sample(B)
n=100
Get <- function(a){
sum(sample(A,a))*Wa + sum(sample(B,n-a))*Wb
}
aa=0:n
k=10^4
MAX=replicate(k,which.max(sapply(aa,Get)))
summary(MAX)
hist(MAX,freq=FALSE)
lines(density(MAX),lwd=2)
132:卵の名無しさん
17/10/08 06:36:02.71 HZvOcnfD.net
>>128
パチンコのスレにこういう記載があった。
>バカは何故かパチンコで勝てると思い込んでいる
>バカは本来なら勝てるのに遠隔や不正で負けていると思い込んでいる
>バカは10分の1で10円が当たるクジより1万分の1で8000円が当たるクジの方が儲かると思ってる
>算数すらできないバカが必死に守って支えてきたのがパチンコ業界
これを読んでこんな問題を考えてみた。
宝くじAは1000本中100本が当たりで当たりは 100万円の賞金、
宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、
どちらも売り出し価格は同じなので100本買うことにする。
どちらを何本買うときに賞金の期待値が最大になるか、シミュレーションしてみよ。
133:卵の名無しさん
17/10/08 12:43:27.27 jiwdTQ4W.net
MAP <- function(x) { # モード値を計算
dens <- density(x)
mode_i <- which.max(dens$y)
mode_x <- dens$x[mode_i]
mode_y <- dens$y[mode_i]
c(x=mode_x, y=mode_y)
}
Na=1000 ; Nb=1000
pa=0.1 ; pb=0.01
Wa=100 ; Wb=1000
n=100
A0=c(rep(1,Na*pa),rep(0,Na-Na*pa))
B0=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb))
Get <- function(a){ #Aをa本買ったときの賞金合計
sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb
}
takarakuji <- function(k=10^3,Histogram=FALSE,Print=TRUE,...){
MAX=replicate(k,which.max(sapply(aa,Get))-1)
if(Histogram) {hist(MAX,freq=FALSE,...) ; lines(density(MAX),lwd=2)}
if(Print)print(round(c(mean=mean(MAX),sd=sd(MAX),mode=MAP(MAX)[1]),2))
invisible(c(mean=mean(MAX),sd=sd(MAX),mode=MAP(MAX)[1]))
}
takarakuji(10^5,Histogram = TRUE,col='lightblue')
134:卵の名無しさん
17/10/08 18:38:07.36 jiwdTQ4W.net
Takarakuji <- function(k=10^2,.Na=1000,.Nb=1000,.pa=0.1,.pb=0.01,
.n=100,.Wa=100,.Wb=1000,Histogram=TRUE,Print=TRUE,...){
Na=.Na ; Nb=.Nb
pa=.pa ; pb=.pb
Wa=.Wa ; Wb=.Wb
n=.n
aa=0:n
A0=c(rep(1,Na*pa),rep(0,Na-Na*pa))
B0=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb))
Get <- function(a)sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb
MAX=replicate(k,which.max(sapply(aa,Get))-1)
tbl_MAX=table(MAX)
MODE=tbl_MAX[which.max(tbl_MAX)]
mode=as.numeric(names(MODE))
if(Histogram) {
hist(MAX,freq=FALSE,main=paste('Na=',Na,'Nb=',Nb),...)
lines(density(MAX),lwd=1)
}
if(Print){
print(round(c(mean=mean(MAX),sd=sd(MAX),mode=mode),2))
}
invisible(c(mean=mean(MAX),sd=sd(MAX),mode=mode))
}
dev.off()
par(mfrow=c(2,2))
Takarakuji(.Na=10^3,.Nb=10^3,col='lightblue')
Takarakuji(.Na=10^5,.Nb=10^5,col='wheat')
Takarakuji(.Na=10^7,.Nb=10^7,col='maroon')
Takarakuji(.Na=10^8,.Nb=10^8,col='gray')
135:卵の名無しさん
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)