「安倍が死ねば、日本は幸せになる」という命題があるとき、安倍が死んでないならこの命題は必ず真 [無断転載禁止]©2ch.net

1132人目の素数さん2017/03/29(水) 06:03:13.82ID:LYDp2KuN
これマジ?

170◆2VB8wsVUoo 2017/09/08(金) 17:11:13.71ID:6ibQhXIy

171◆2VB8wsVUoo 2017/09/08(金) 17:11:30.22ID:6ibQhXIy

172◆2VB8wsVUoo 2017/09/08(金) 17:11:48.12ID:6ibQhXIy

173◆2VB8wsVUoo 2017/09/08(金) 17:12:05.47ID:6ibQhXIy

174◆2VB8wsVUoo 2017/09/08(金) 17:12:23.02ID:6ibQhXIy

175◆2VB8wsVUoo 2017/09/08(金) 17:12:41.22ID:6ibQhXIy

176◆2VB8wsVUoo 2017/09/08(金) 17:12:59.84ID:6ibQhXIy

177◆2VB8wsVUoo 2017/09/08(金) 17:13:19.29ID:6ibQhXIy

178◆2VB8wsVUoo 2017/09/08(金) 17:13:38.03ID:6ibQhXIy

179132人目の素数さん2017/09/09(土) 12:20:51.73ID:G2DuD1v6
耳栓をしたら世界が変わってワロタ

180◆2VB8wsVUoo 2017/10/24(火) 08:24:40.08ID:n4VIAFJn

181◆2VB8wsVUoo 2017/10/24(火) 08:24:58.28ID:n4VIAFJn

182◆2VB8wsVUoo 2017/10/24(火) 08:25:15.50ID:n4VIAFJn

183◆2VB8wsVUoo 2017/10/24(火) 08:25:31.64ID:n4VIAFJn

184◆2VB8wsVUoo 2017/10/24(火) 08:25:49.57ID:n4VIAFJn

185◆2VB8wsVUoo 2017/10/24(火) 08:26:07.32ID:n4VIAFJn

186◆2VB8wsVUoo 2017/10/24(火) 08:26:26.09ID:n4VIAFJn

187◆2VB8wsVUoo 2017/10/24(火) 08:26:50.38ID:n4VIAFJn

188◆2VB8wsVUoo 2017/10/24(火) 08:27:09.09ID:n4VIAFJn

189◆2VB8wsVUoo 2017/10/24(火) 08:27:29.85ID:n4VIAFJn

190132人目の素数さん2018/04/07(土) 07:01:10.55ID:eG1yamTB
jags4prop <- function(r1,r2,n1,n2,
ROPEdiff=c(-0.025,0.025),ROPErate=c(0.80,1.25),
NNT=FALSE){
library(rjags)
y=c(rep(1,r1),rep(0,n1-r1),rep(1,r2),rep(0,n2-r2))
s=as.numeric(factor(c(rep('D',n1),rep('U',n2))))
a=1 ; b=1 # JAGS prior : beta(a,b)
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
# 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," )
}
}
")
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
mcmcMat=as.matrix(codaSamples)
risk1=mcmcMat[,1]
risk2=mcmcMat[,2]
x11()
par(mfrow=c(2,2))
BEST::plotPost(risk1,col='gray')
if(NNT){
BEST::plotPost(abs(1/(risk1-risk2)), xlab ='NNT(NNH)',col='pink')}
else{
BEST::plotPost(risk1-risk2,compVal=0,ROPE=ROPEdiff,cex=1,col='pink')
}
BEST::plotPost(risk2,col='gray')
BEST::plotPost(risk1/risk2,compVal=1,ROPE=ROPErate,col='skyblue')
invisible(codaSamples)
}
jags4prop(r1=1,n1=1,r2=0,n2=1,ROPEdiff=NULL,NNT=TRUE)
jags4prop(r1=1,n1=1,r2=0,n2=1,ROPEdiff=NULL)

191132人目の素数さん2018/04/11(水) 11:58:17.41ID:LX0qXOLG
f<- function(d=1,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f(5/10)
power.t.test(delta=0.5,power=0.80)$n
# power 50%->80%
dd=seq(0.01,0.99,by=0.01)
g<- function(x) f(x,.50)/f(x,.80)
summary (sapply(dd,g))
(2*(abs(1.96+abs(qnorm(1-0.80))))^2)/(2*(abs(1.96+abs(qnorm(1-0.50))))^2)

(1.96+0.84)^2/(1.96+0)^2

dd=seq(0.01,0.99,by=0.01)
h<-function (x,pow1=0.80,pow2=0.50){
power.t.test(d=x, power=pow1)$n/power.t.test(d=x, power=pow2)$n
}
res=sapply (dd,h)
plot(dd,res)
summary (res)

# power 80%->90%
f(beta=0.90)
f(beta=0.80)
(2*(abs(1.96+abs(qnorm(1-0.90))))^2)/(2*(abs(1.96+abs(qnorm(1-0.80))))^2)
res1=sapply (dd, function (x) h(x,0.90,0.80))
plot (dd,res1)

# power 80%->99%
f(beta=0.99)
f(beta=0.80)
(2*(abs(1.96+abs(qnorm(1-0.99))))^2)/(2*(abs(1.96+abs(qnorm(1-0.80))))^2)
res2=sapply (dd, function (x) h(x,0.99,0.80))
plot (dd,res2)

##
f<- function(d=1,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f()
sample.size <- function (d=1,beta=0.80,alpha=0.05){
f<- function(n) 2*((abs(qt(1-alpha/2,n-1))+abs(qt(1-beta,n-1)))/d)^2 - n
uniroot (f, c(3,10000))$root
}
sample.size()
power.t.test(d=1, power=0.80)$n 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)

192132人目の素数さん2018/04/19(木) 23:50:23.50ID:GVMUXyX9
慶應義塾学医学部 2051 70
順天堂大学医学部 2090 66
昭和大学医学部 2200 65
東京慈恵会医科大学医学部 2250 69
自治医科大学医学部 2260 66
産業医科大学医学部 2749 65
日本医科大学医学部 2813 66
東京医科大学医学部 2995 64
関西医科大学医学部 3014 65
大阪医科大学医学部 3141 66
東邦大学医学部 3180 62
久留米大学医学部 3237 61
東京女子医科大学医学部 3284 60
日本大学医学部 3310 61
岩手医科大学医学部 3400 60
聖マリアンナ医科大学医学部 3440 58
近畿大学医学部 3580 62
藤田保健衛生大学医学部 3650 60
獨協医科大学医学部 3660 59
帝京大学医学部 3750 59
杏林大学医学部 3755 61
東海大学医学部 3760 61
福岡大学医学部 3770 60
愛知医科大学医学部 3800 59
埼玉医科大学医学部 3800 57
兵庫医科大学医学部 3880 60
北里大学医学部 3890 60
金沢医科大学医学部 3950 58
川崎医科大学医学部 4565 57

dat=read.table("clipboard")
colnames(dat)=c("school","tuition","hensa")
head(dat)

xi=dat$hensa
yi=dat$tuition
re=lm(yi~xi) ; re
beta1.hat=re$coef[1]
beta10=0
n=length(xi)
Sxx=var(xi)*(n-1)
Syy=var(yi)*(n-1)
Sxy=cov(xi,yi)*(n-1)
beta1.hat=Sxy/Sxx
Se=Syy - Sxy^2/Sxx
phie=n-2
Ve=Se/phie
t0=(beta1.hat-beta10)/sqrt(Ve/Sxx) ; t0
curve(dt(x,n-2),-5,5,bty='l')
pt(t0,n-2)

> t0=(beta1.hat-beta10)/sqrt(Ve/Sxx) ; t0
[1] -10.17626
> pt(t0,n-2)
[1] 4.873385e-11

193132人目の素数さん2018/04/27(金) 20:41:28.71ID:uC0CNGcQ
# Lindner Center data on 996 PCI patients analyzed by Kereiakes et al. (2000)

library(MatchLinReg)
data("lindner")
head(lindner)
attach(lindner)
(tbl=table(stent,sixMonthSurvive))
Epi::twoby2(tbl)
formula.stent = stent ~ abcix + height + female + diabetic + acutemi +
ejecfrac + ves1proc - 1
formula.6month = sixMonthSurvive ~ stent+ abcix + height + female + diabetic + acutemi +
ejecfrac + ves1proc - 1
rms::lrm( stent ~ abcix + height + female + diabetic + acutemi +
ejecfrac + ves1proc-1, data=lindner)$stat['C']
ps=re.glm$fitted.values
Y=lindner$sixMonthSurvive
Tr=lindner$stent

# IPWE
weighted.mean(Y,Tr/ps) - weighted.mean(Y,(1-Tr)/(1-ps))

## (Doubly Robust: DR)
dre <- function(data, target, treat, ps, formula) {
n <- nrow(data)
y <- data[target]
z <- data[treat]
data1 <- data[data[treat]==1,]
data0 <- data[data[treat]==0,]
model1 <- lm(formula=formula, data=data1)
model0 <- lm(formula=formula, data=data0)
fitted1 <- predict(model1, data)
fitted0 <- predict(model0, data)
dre1 <- (1/n)*sum(y+((z-ps)/ps)*(y-fitted1))
dre0 <- (1/n)*sum(((1-z)*y)/(1-ps)+(1-(1-z)/(1-ps))*fitted0)
c(treated=dre1,control=dre0,ATE=dre1-dre0)
}

dre(data=lindner,target='sixMonthSurvive',treat='stent', ps=ps, formula=sixMonthSurvive ~ abcix + height
+ female + diabetic + acutemi +ejecfrac + ves1proc -1)

194132人目の素数さん2018/05/03(木) 11:38:00.90ID:Jusm5aoQ
# http://www1.udel.edu/educ/gottfredson/reprints/2009semen.pdf
cor2ci<- function(cor,n,conf.level=0.95){
lower=tanh(atanh(cor) + qnorm((1-conf.level)/2)/sqrt(n-3))
upper=tanh(atanh(cor) - qnorm((1-conf.level)/2)/sqrt(n-3))
data.frame(lower,cor,upper)
}

cor2p <- function(cor,n,alternative='two.sided'){
r=cor
t.stat=r*sqrt(n-2)/sqrt(1-r^2)
if(alternative=='two.sided'){
p.value=pt(-abs(t.stat),n-2)*2
}else{
p.value=pt(-abs(t.stat),n-2)
}
return(p.value)
}
cor2p(0.14,425)
rs=c(conc=0.15,count=0.19,Motility=0.14)
n=425
cor2ci(rs,n)
sapply(rs,function(cor)cor2p(cor,n=425))
#
r2p <- function(r,n){
2*pt(-abs(r)*sqrt(n-2)/sqrt(1-r^2),n-2)
}
# 0.002=r2p(0.15,n)
curve(r2p(rs[1],x),300,600)
abline(h=0.002,lty=3)
uniroot(function(n,p0=0.0015)r2p(rs[1],n)-p0,c(100,1000))$root # 445
uniroot(function(n,p0=0.0020)r2p(rs[1],n)-p0,c(100,1000))$root # 422
uniroot(function(n,p0=0.0025)r2p(rs[1],n)-p0,c(100,1000))$root # 404
r2p(rs[1],425)

# 0.002=r2p(0.14,n)
curve(r2p(rs[3],x),300,600)
abline(h=0.002,lty=3)
uniroot(function(n,p0=0.0015)r2p(rs[3],n)-p0,c(100,1000))$root # 511
uniroot(function(n,p0=0.0020)r2p(rs[3],n)-p0,c(100,1000))$root # 484
uniroot(function(n,p0=0.0025)r2p(rs[3],n)-p0,c(100,1000))$root # 464
r2p(rs[3],425)
r2p(rs[3],c(445,422,404))

r2p(0.19,425)

195132人目の素数さん2018/05/04(金) 20:19:11.93ID:Za4pacSv
noise2sig <− function(seed,sequential=FALSE){
set.seed(seed)
Noise=matrix(rnorm(100*51),ncol=51,nrow=100)
# first pass
lm1=lm(Noise[,51] 〜 Noise[,−51]) # Noise[,51] : 目的変数Y
cat(V\nR2 (1st) = V,summary(lm1)$r.squared)
cat(V\np.value (1st) = V,anova(lm1)$Pr[1])
coefs1=summary(lm1)$coef[,VPr(>|t|)V]
length(coefs1)
coefs1[1]
coeffs1=coefs1[−1]
cat(V\ncoef < 0.25 =V,sum(coeffs1<0.25))
cat(V\ncoef < 0.05 =V,sum(coeffs1<0.05))
indx25=which(coeffs1<0.25)

# second pass
lm2=lm(Noise[,51] 〜 Noise[,indx25])
cat(V\n\nR2 (2nd) = V,summary(lm2)$r.squared)
cat(V\np.value(2nd) = V,anova(lm2)$Pr[1])
coefs2=summary(lm2)$coef[,VPr(>|t|)V]
length(coefs2)
coefs2[1]
coeffs2=coefs2[−1]
cat(V\ncoef < 0.25 (2nd) =V,sum(coeffs2<0.25))
cat(V\ncoef < 0.05 (2nd) =V,sum(coeffs2<0.05))
cat(V\n\nV)
if(sequential){
  cat(VHit Return Key in console windowV)
  no_save <− scan(n=1, what=character(), quiet=TRUE)
 }



noise2sig(1)
for(i in 2:5) noise2sig(i,seq=TRUE)

196132人目の素数さん2018/05/08(火) 13:22:38.87ID:Gebq//ip
> UG
1年度 2年度 3年度 4年度 5年度 6年度
学生1 30 45 42 31 47 35
学生2 42 18 46 42 42 52
学生3 51 39 23 45 44 41
学生4 28 56 26 49 41 33
学生5 51 48 31 48 42 37
学生6 41 58 27 37 38 30
学生7 20 31 42 16 31 29
学生8 45 42 54 55 46 50
学生9 61 35 34 38 47 34
学生10 44 32 23 30 29 56
> mUG=reshape2::melt(UG) ; colnames(mUG)=c('学生','年度','成績')

> Tk=TukeyHSD(aov(成績~学生,mUG))

> p.Tk=Tk[[1]][,'p adj']
> n9=nrow(UG)-1
> mat.Tk=matrix(rep(NA,n9*n9),n9)
> for(k in 1:9) mat.Tk[,k]=c(rep(NA,k-1), p.Tk[((k-1)*(n9+1-k/2)+1):((k-1)*(n9+1-k/2)+1+n9-k)])
> colnames(mat.Tk)=paste0('学生',1:9)
> rownames(mat.Tk)=paste0('学生',2:10)
> d=round(mat.Tk,3)
> d[is.na(d)]='-'
> print(d,quote=FALSE)
学生1 学生2 学生3 学生4 学生5 学生6 学生7 学生8 学生9
学生2 1 - - - - - - - -
学生3 1 1 - - - - - - -
学生4 1 1 1 - - - - - -
学生5 0.998 1 1 0.999 - - - - -
学生6 1 1 1 1 0.999 - - - -
学生7 0.742 0.513 0.493 0.687 0.257 0.724 - - -
学生8 0.724 0.9 0.911 0.776 0.989 0.742 0.023 - -
学生9 1 1 1 1 1 1 0.383 0.958 -
学生10 1 0.998 0.997 1 0.958 1 0.945 0.418 0.989
>

197132人目の素数さん2018/07/22(日) 11:22:57.29ID:Ott8rTSz
> REDWHITECOLOR(perm[indx,])
RED = 5 8 7 WHITE = 3 9 6 1 8 COLOR = 4 0 2 0 5

198132人目の素数さん2018/07/22(日) 11:43:12.65ID:Ott8rTSz
# ド底辺+シリツ=裏口馬鹿

x=LETTERS[1:10]
uraguchi <- function(A,B,C,D,E,F,G,H,I,J){
doteihen=10^(2:0)
siritsu=10^(2:0)
uraguchibaka=10^(3:0)
sum(doteihen*c(A,B,C))+sum(siritsu*c(D,E,F))-sum(uraguchibaka*c(G,H,I,J))
}
x=unique(x)
URAGUCHI <- function(x){
A=x[1]
B=x[2]
C=x[3]
D=x[4]
E=x[5]
F=x[6]
G=x[7]
H=x[8]
I=x[9]
J=x[10]

cat(paste('ド底辺 = ',A,B,C,
' シリツ = ',D,E,F,
' 裏口馬鹿 = ',G,H,I,J),'\n')
}
URAGUCHI(unique(x))

library(gtools)
perm=permutations(n=10,r=10,v=0:9)
perm=perm[perm[,1]!=0&perm[,4]!=0&perm[,7]!=0,] # A!=0,D!=0,G!=0
head(perm) ; tail(perm)
n=nrow(perm)
re=numeric(n)
for(i in 1:n){
re[i]=uraguchi(perm[i,1],perm[i,2],perm[i,3],perm[i,4],perm[i,5],perm[i,6],perm[i,7],perm[i,8],perm[i,9],perm[i,10])
}
hist(re)
indx=which(re==0)
(n.indx=length(indx))
for(i in 1:n.indx) URAGUCHI(perm[indx[i],])

199132人目の素数さん2018/07/22(日) 15:29:39.25ID:Ott8rTSz
/*
  底辺私立
 ×   5
------------
  裏口馬鹿
各々の漢字は1〜8までの数字を表します。
同じ数字は使われません。
底辺私立・裏口馬鹿に数字をいれて筆算を2通り完成させてください。
*/

#include<stdio.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;
}
int unique(int num[]){
int i,j,n=8;
qsort(num,n,sizeof(int),compare_int);
for(i=0;i<n;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;
}
main(){
int A,B,C,D,E,F,G,H;
for(A = 1; A < 9; A++){
for(B = 1; B < 9; B++){
for(C = 1; C < 9; C++){
for(D = 1; D < 9; D++){
for(E = 1; E < 9; E++){
for(F = 1; F < 9; F++){
for(G = 1; G < 9; G++){
for(H = 1; H < 9; H++){
if((A*1000+B*100+C*10+D)*5==E*1000+F*100+G*10+H){
int num[]={A,B,C,D,E,F,G,H};
if(unique(num)==1){
printf("底辺私立 = %2d%2d%2d%2d 裏口馬鹿 = %2d%2d%2d%2d\n", A,B,C,E,F,G,H);
}}
}}} }}} }}
}

200132人目の素数さん2018/07/31(火) 12:09:36.10ID:AXaAgLMZ
attack <- function(x){ # TRUE if on attack
n=length(x)
for(i in 1:(n-1)){
for(j in 1:(n-i)){
if(x[i+j]==x[i]+j | x[i+j]==x[i]-j) return(TRUE)
}
}
return(FALSE)
}

N_Queen<-function(n){
perm=gtools::permutations(n,n)
ret=apply(perm,1,attack)
perm[which(ret==FALSE),]
}

N_Queen(8)

201132人目の素数さん2018/07/31(火) 12:11:17.41ID:AXaAgLMZ
rot <- function(row,col,n=8){ # position after 90 degree rotation
c(col,n-row+1)
}
rota <- function(x){ # 90 degree rotation
n=length(x)
re=numeric(n)
for(i in 1:n){
tmp=rot(i,x[i],n)
re[tmp[1]]=tmp[2]
}
re
}

# ratete 90 degree four times
rotate <- function(x){
r1=rota(x)
r2=rota(r1)
r3=rota(r2)
rx=rbind(x,r1,r2,r3)
rownames(rx)=NULL
return(rx)
}

b2v <- function(x){
ret=unique.matrix(rbind(x,rotate(x),rotate(rev(x))))
rownames(ret)=NULL
ret
}

is.var <- function(x,y){ # x: base, y : variation of x ?
v=b2v(x)
nv=nrow(v)
for(i in 1:nv){
if(all(v[i,]==y)) return(TRUE)
}
return(FALSE)
}

uniQ <- function(M){ # check firt row vs the others
n=nrow(M)
uniq=numeric(n)
for(i in 1:n){
uniq[i]=!is.var(M[1,],M[i,]) # is no variant?
}
M=rbind(M[which(uniq==1),],M[1,]) # append 1st row at last
return(M) # not the final answer
}
M=N_Queen(9)
while(nrow(M)!=nrow(uniQ(M))){
M=uniQ(M)
}
M

202132人目の素数さん2018/07/31(火) 13:30:25.25ID:AXaAgLMZ
basic_Queen <- function(N){
M=m=N_Queen(N)
while(nrow(M)!=nrow(uniQ(M))) M=uniQ(M)
return(list(variant=m,basic=M))
}

basic_Queen(5)

drawQ <- function(n){#
M=basic_Queen(n)$basic
nr=nrow(M)
for(i in 1:nr){
cat('(',i,')\n')
for(j in 1:n){
cat(rep('_',M[i,j]-1),'●',rep('_',n-M[i,j]),'\n')
}
cat('\n')
}
}
drawQ(5)
drawQ(6)
drawQ(8)

203132人目の素数さん2018/08/07(火) 11:36:28.79ID:JOZRDgmj
/* ウカウカ+ウサギ+ノソノソ+カメ+キヨウソウ=ウサギトカメ */
#include<stdio.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;
}
int unique(int num[]){
int i,j,n=10;
qsort(num,n,sizeof(int),compare_int);
for(i=0;i<n;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;
}
main(){
int u,ka,sa,gi,no,so,me,ki,yo,to;
for(u = 1; u < 10; u++){
for(ka = 1; ka < 10; ka++){
for(sa = 0; sa < 10; sa++){
for(gi = 0; gi < 10; gi++){
for(no = 1; no < 10; no++){
for(so = 0; so < 10; so++){
for(me = 0; me < 10; me++){
for(ki = 1; ki < 10; ki++){
for(yo = 0; yo < 10; yo++){
for(to = 0; to < 10; to++){
/* ウカウカ+ウサギ+ノソノソ+カメ+キヨウソウ=ウサギトカメ */
if(u*1000+ka*100+u*10+ka +u*100+sa*10+gi +no*1000+so*100+no*10+so +ka*10+me +ki*10000+yo*1000+u*100+so*10+u == u*100000+sa*10000+gi*1000+to*100+ka*10+me){
int num[]={u,ka,sa,gi,no,so,me,ki,yo,to};
if(unique(num)==1){
printf("%d%d%d%d + %d%d%d + %d%d%d%d + %d%d + %d%d%d%d%d = %d%d%d%d%d%d\n",u,ka,u,ka,u,sa,gi,no,so,no,so,ka,me,ki,yo,u,so,u,u,sa,gi,to,ka,me);
}}
}}}}}}}}}}
} 👀
Rock54: Caution(BBR-MD5:1341adc37120578f18dba9451e6c8c3b)

204132人目の素数さん2018/08/10(金) 12:25:59.56ID:lOg+llmH
# 7個の a と3個の b を一列に並べてできる順列のうち
# 次の簡約律のもとで文字を消していくと最終的に何も残らなくなる順列は何通りありますか。
# ・aa が現れると消える。
# ・bb が現れると消える。
# ・ababab が現れると消える。
# ・bababa が現れると消える。

# all permutation
reperm <- function(n, r, v) {
# if(n!=length(v) | n < r) return('ERROR')
if(r==1) return(matrix(v))
else {
re <- NULL
for(i in 1:n) re=rbind(re, cbind(v[i], reperm(n-1,r-1,v[-i])))
return(unique(re))
}
}
perm=reperm(n=10,r=10,v=c(rep(0,7),rep(1,3)))

# ababab|bababa
rm_ababab <- function(v){
if(length(v)<6) return(v)
for(i in 1:5)
if(all(v[i:(i+5)]==c(0,1,0,1,0,1))|all(v[i:(i+5)]==c(1,0,1,0,1,0))){
return(v[-i:-(i+5)])
}
return(v)
}

# aa|bb
rm_aa <- function(v){ # remove 1st two aa or bb
if(length(v)<2) return(v)
n=length(v)
for(i in 1:(n-1))
if(all(v[i:(i+1)]==c(0,0))|all(v[i:(i+1)]==c(1,1))){
return(v[-i:-(i+1)])
}
return(v)
}

rm_aaR <- function(v){ # recursively remove aa or bb while checking abababa or bababa
if(length(v)==6){
if(all(v==c(1,0,1,0,1,0))|all(v==c(0,1,0,1,0,1))) return(c(0,0))
}
if(length(v)==2) return(v)
else{
v1=rm_aa(v)
Recall(v1)
}
}

res=apply(perm,1,rm_aaR)
res=t(res)
sum(res[,1]==0 & res[,2]==0)
idx=which(res[,1]==0 & res[,2]==0)
perm[idx,]

205132人目の素数さん2018/08/10(金) 15:19:15.44ID:Hlm8Oe3x
combnを使って改訂
# all permutations
ones=t(combn(10,3))
indx2ten <- function(x){
y=numeric(10)
y[x[1]]=1
y[x[2]]=1
y[x[3]]=1
return(y)
}
perm=NULL
for(i in 1:nrow(ones)) perm=rbind(perm,indx2ten(ones[i,]))
# ababab|bababa
rm_ababab <- function(v){
if(length(v)<6) return(v)
for(i in 1:5)
if(all(v[i:(i+5)]==c(0,1,0,1,0,1))|all(v[i:(i+5)]==c(1,0,1,0,1,0))){
return(v[-i:-(i+5)])
}
return(v)
}

# aa|bb
rm_aa <- function(v){ # remove 1st two aa or bb
if(length(v)<2) return(v)
n=length(v)
for(i in 1:(n-1))
if(all(v[i:(i+1)]==c(0,0))|all(v[i:(i+1)]==c(1,1))){
return(v[-i:-(i+1)])
}
return(v)
}
rm_aaR <- function(v){ # recursively remove aa or bb while checking abababa or bababa
if(length(v)==6){
if(all(v==c(1,0,1,0,1,0))|all(v==c(0,1,0,1,0,1))) return(c(0,0))
}
if(length(v)==2) return(v)
else{
v1=rm_aa(v)
Recall(v1)
}
}
res=apply(perm,1,rm_aaR)
res=t(res)
sum(res[,1]==0 & res[,2]==0)
idx=which(res[,1]==0 & res[,2]==0)
Ans=perm[idx,]
AB=Ans+1
rownames(AB)=NULL
.ab=c('a','b')
indx2char <- function(x,ab=.ab){
n=length(x)
re=NULL
for(i in 1:n) re[i]=ab[x[i]]
re
}
print(t(apply(AB,1,indx2char)),quote = FALSE)

206132人目の素数さん2018/09/04(火) 13:44:16.80ID:mJ0F6GVN
A:安倍が死ぬ
B:日本は幸せ

A∧(A→B))

A∧(¬A∨B)

¬(A∧¬A)

B

207132人目の素数さん2018/09/11(火) 13:59:25.72ID:KvhdapkQ
Gacha.fm <- function(p,write=FALSE){
n=length(p)
par=letters[1:n]
fm <- function(v){
nv=length(v)
re=character(nv)
for(j in 1:nv) re[j]=par[v[j]]
s=paste(re,collapse='+')
if(nv==1) paste0('1/',s)
else paste0('1/(',s,')')
}
fm1 <- function(mat){
paste(apply(mat,2,fm),collapse='+')
}
re=list()
for(i in 1:n) re[[i]]=fm1(combn(n,i))
re1=re[[1]]
for(i in 2:(n-1)){
re1=c(re1,ifelse(i%%2,' + ',' - '),'{',re[[i]],'}')
}
output=paste(paste(re1,collapse=""),ifelse(n%%2,'+','-'), re[[n]])
cat(output,'\n')
if(write) write(output,'output.txt')
invisible(output)
}

208132人目の素数さん2018/10/26(金) 16:44:10.02ID:Jik/lAlw
jags4prop <- function(r1,r2,n1,n2,
ROPEdiff=c(-0.025,0.025),ROPErate=c(0.80,1.25),
NNT=FALSE){
library(rjags)
y=c(rep(1,r1),rep(0,n1-r1),rep(1,r2),rep(0,n2-r2))
s=as.numeric(factor(c(rep('D',n1),rep('U',n2))))
a=1 ; b=1 # JAGS prior : beta(a,b)
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
# 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," )
}
}
")
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
mcmcMat=as.matrix(codaSamples)
rate1=mcmcMat[,1]
rate2=mcmcMat[,2]
x11()
par(mfrow=c(2,2))
BEST::plotPost(rate1,col='gray')
if(NNT){
BEST::plotPost(abs(1/(rate1-rate2)), xlab ='NNT(NNH)',col='pink')}
else{
BEST::plotPost(rate1-rate2,compVal=0,ROPE=ROPEdiff,cex=1,col='pink')
}
BEST::plotPost(rate2,col='gray')
BEST::plotPost(rate1/rate2,compVal=1,ROPE=ROPErate,col='skyblue')
invisible(codaSamples)
}

jags4prop(50,19,133,113)

209132人目の素数さん2018/10/30(火) 00:30:28.16ID:TZqGbv4d
--コンパイルしてコマンドラインから実行できるように改変(但し、エラー処理皆無)

import System.Environment
import Data.List
import Data.List.Split

combinations :: Int -> [a] -> [[a]]
combinations 0 _ = [ [] ]
combinations n xs = [ y:ys | y:xs' <- tails xs, ys <- combinations (n-1) xs']

main = do
argList <- getArgs -- m : 縦マス(短軸) n : 横マス(長軸) k : 宝の数
let m = read (argList !! 0)
n = read (argList !! 1)
k = read (argList !! 2)
q = [0..m*n-1]
matQ = chunksOf n q
matP = transpose matQ --行列を転置して
p = concat matP -- 配列に変換
treasure = combinations k q -- 宝の組み合わせ
ip y = minimum $ map(\x -> elemIndices x p!!0) y -- 宝の、配列pでのindex列を求めて最小値を返す
iq y = minimum $ map(\x -> elemIndices x q!!0) y
idxp = map ip treasure -- 宝の組み合せで実行して
idxq = map iq treasure
p_q = zipWith (-) idxp idxq -- 差をとって大小判別
p1st = length $ filter (<0) p_q -- 短軸方向探索pが先に宝をみつける
q1st = length $ filter (>0) p_q
draw = length $ filter (==0) p_q
putStrLn $ "p1st = " ++ show p1st ++ ", q1st = " ++ show q1st ++ ", draw = " ++ show draw

210132人目の素数さん2018/11/09(金) 20:08:31.89ID:5AnUTlVm
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95){
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)
invisible(ICDF)
}

# N個のクジでr個めで初めてあたった時のN個内の当たり数の推測
Atari <- function(N,r){
pmf <- function(x) ifelse(x>N-r+1,0,(1-x/N)^(r-1)*x/N) # dnbinom(r-1,1,x/N) ; dgeom(r-1,x/N)
# curve((1-x/N)^(r-1)*x/N,0,N)
AUC=integrate(pmf,0,N)$value
pdf <- function(x) pmf(x)/AUC
mode=optimise(pdf,c(0,N),maximum=TRUE)$maximum
mean=integrate(function(x)x*pdf(x),0,N)$value
cdf <- function(x) integrate(pdf,0,x)$value
median=uniroot(function(x)cdf(x)-0.5,c(0,N))$root
print(c(mode=mode,median=median,mean=mean))
pdf2hdi(pdf,0,N,cred=0.95)
}
Atari(100,3)

# 全体N個中当たりS個、1個ずつ籤を引いて当たったらやめる.
# r個めが初めて当たりであったときSの信頼区間を推定するシミュレーション。

atari <- function(N,r,k=1e3){ # k: simlation times
f <- function(S,n=N){which.max(sample(c(rep(1,S),rep(0,n-S))))}
vf=Vectorize(f)
sim=replicate(k,vf(1:(N-r)))
s=which(sim==r)%%(N-r)
s[which(s==0)]=N-r
hist(s,freq=T,col='skyblue')
print(quantile(s,c(.025,.05,.50,.95,.975)))
print(HDInterval::hdi(s))
}
atari(100,3,1e3)

211132人目の素数さん2018/11/11(日) 20:15:20.63ID:ODPKEEGK
from fractions import Fraction

def seq_dice(N,k,p=1/2):
P=list()
for n in range(k-1):
P.append(0)
P.append(p**k)
P.append(p**k + (1-p)*p**k)
for n in range (k,N):
P.append(P[n]+(1-P[n-k])*p**(k+1))
return(P[N])

def seq_diceJ(N,k,p=1/2):
return(seq_dice(N,k,p) - seq_dice(N,k+1,p))


def dice(N,k,p=1/2):
print("Over " + str(k))
print(Fraction(seq_dice(N,k,p)))
print(" = " + str(seq_dice(N,k,p)))
print("Just " + str(k))
print(Fraction(seq_diceJ(N,k,p)))
print(" = " + str(seq_diceJ(N,k,p)))


dice(100000000,25)

212132人目の素数さん2018/11/11(日) 20:27:20.16ID:ODPKEEGK
Over 25
6977459029519597/9007199254740992
= 0.7746535667951356
Just 25
2246038053550679/9007199254740992
= 0.24936031612362342

213132人目の素数さん2018/11/11(日) 20:27:38.50ID:ODPKEEGK
seqNp <- function(N=100,K=5,p=0.5){
if(p==0) return(0)
if(N==K) return(p^K)
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q
for(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}
P0=numeric(N)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0
for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)
ret[N]
}

seqNpJ <- function(N,K,p) seqNp(N,K,p)-seqNp(N,K+1,p)
seqNpJ(100,5,0.5)

214132人目の素数さん2018/11/11(日) 20:28:04.68ID:ODPKEEGK
seq2pCI <- function(N,K,alpha=0.05,Print=T){
vp=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
if(Print){curve(vp(x),lwd=2,bty='l',xlab='Pr[head]',ylab=paste('Pr[max',K,'-head repetition]'))
abline(h=alpha,lty=3)}
peak=optimize(vp,c(0,1),maximum=TRUE)$maximum
mean=integrate(function(x)x*vp(x),0,1)$value/integrate(function(x)vp(x),0,1)$value
lwr=uniroot(function(x,u0=alpha) vp(x)-u0,c(0.01,peak))$root
upr=uniroot(function(x,u0=alpha) vp(x)-u0,c(peak,0.99))$root
c(lower=lwr,mean=mean,mode=peak,upper=upr)
}
seq2pCI(100,5,0.05,T) # pmfからの計算で正しいか?

vs=Vectorize(function(K)seq2pCI(N=100,K,alpha=0.05,Print=F))
x=2:23
y=vs(x)
head(y)
y=y*100
plot(x,y['mean',],bty='l',pch=19,ylim=c(0,100),
xlab="最大連続数",ylab="推定裏口入学者数")
points(2:23,y['mode',],bty='l')
segments(x,y['lower',],x,y['upper',])
legend('right',bty='n',pch=c(19,1),legend=c("期待値","最頻値"))

215132人目の素数さん2018/11/11(日) 20:29:05.25ID:ODPKEEGK
# pdfからcdfの逆関数を作ってhdiを表示させて逆関数を返す
# 0,1での演算を回避 ∫[1/nxxx,1-1/nxx]dxで計算
pdf2HDI <- function(pdf,xMIN=0,xMAX=1,cred=0.95,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)
c(hdi[1],hdi[2])
}

# N試行で最大K回連続成功→成功確率pの期待値、最頻値と95%HDI
# max K out of N-trial to probability & CI
mKoN2pCI <- function(N=100 , K=4 , conf.level=0.95){
pmf=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
mode=optimize(pmf,c(0,1),maximum=TRUE)$maximum
auc=integrate(pmf,0,1)$value
pdf=function(x) pmf(x)/auc
mean=integrate(function(x)x*pdf(x),0,1)$value
curve(pdf(x),lwd=3,bty='l',xlab='Pr[head]',ylab='density')
lu=pdf2HDI(pdf,cred=conf.level)
curve(pdf(x),lu[1],lu[2],type='h',col='lightblue',add=T)
c(lu[1],mean=mean,mode=mode,lu[2])
}
mKoN2pCI(100,4) # one minute plz
mKoN2pCI(100,5)

216132人目の素数さん2018/11/11(日) 20:29:23.48ID:ODPKEEGK
# N試行で最大K回連続成功→成功確率pの期待値、最頻値と95% Quantile
# max K out of N-trial to probability & CI_quantile
mKoN2pCIq <- function(N=100 , K=4 , alpha=0.05){
pmf=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
mode=optimize(pmf,c(0,1),maximum=TRUE)$maximum
auc=integrate(pmf,0,1)$value
pdf=function(x) pmf(x)/auc
curve(pdf(x),bty='l')
mean=integrate(function(x)x*pdf(x),0,1)$value
cdf=function(x) MASS::area(pdf,0,x)
vcdf=Vectorize(cdf)
lwr=uniroot(function(x)vcdf(x)-alpha/2,c(0,mode))$root
upr=uniroot(function(x)vcdf(x)-(1-alpha/2),c(mode,1))$root
c(lower=lwr,mean=mean,mode=mode,upper=upr)
}
mKoN2pCIq(100,4)
mKoN2pCIq(100,5)

## simulation
mKoN2pCIga<-function(N=100,K=4,Print=T){
pmf=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
auc=integrate(pmf,0,1)$value
pdf=function(x) pmf(x)/auc
(mode=optimize(pmf,c(0,1),maximum=TRUE)$maximum)
(mean=integrate(function(x)x*pdf(x),0,1)$value)
(var=integrate(function(x)(x-mean)^2*pdf(x),0,1)$value)
(sd=sqrt(var))
#正規分布で近似してみる
if(Print){
curve(pdf,bty='l',col='red')
curve(dnorm(x,mean,sd),add=T,col='blue')}
c(lower=qnorm(0.025,mean,sd),mode=mode,mean=mean,upper=qnorm(0.975,mean,sd))
}
mKoN2pCIga(100,5)
vmK=Vectorize(function(K) mKoN2pCIga(100,K))
(y=cbind(0,vmK(2:30)))
plot(1:30,y['mean',],pch=19)
points(1:30,y['mode',])

217132人目の素数さん2018/11/11(日) 20:30:07.92ID:ODPKEEGK
# random numbers following PDF by John von Neuman's method
vonNeumann <- function(PDF,xmin=0,xmax=1,N=10000,Print=TRUE,...){
xx=seq(xmin,xmax,length=N+1)
ymax=max(PDF(xx))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
if(Print){
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,col=sample(colors(),1),main='',...)
AUC=integrate(PDF,xmin,xmax)$value
lines(xx,sapply(xx,function(x)PDF(x)/AUC))
}
hdi=HDInterval::hdi(Rand)
print(c(hdi[1],hdi[2]),digits=4)
invisible(Rand)
}

vonNeumann2 <- function(PDF,xmin=0,xmax=1,N=10000,Print=TRUE,...){
xx=seq(xmin,xmax,length=N+1)
xx=xx[-(N+1)]
xx=xx[-1]
ymax=max(PDF(xx))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
if(Print){
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,col=sample(colors(),1),main='',...)
AUC=integrate(PDF,xmin,xmax)$value
lines(xx,sapply(xx,function(x)PDF(x)/AUC))
}
hdi=HDInterval::hdi(Rand)
print(c(hdi[1],hdi[2]),digits=5)
invisible(Rand)
}

218132人目の素数さん2018/11/12(月) 14:19:30.16ID:boi8bvdM
seq_dice <- function(N=100,k=5,p=1/2){
P=numeric(N)
for(n in 1:(k-1)){
P[n]=0
}
P[k]=p^k
P[k+1]=p^k+(1-p)*p^k
for(n in (k+1):(N-1)){
P[n+1] = P[n] + (1-P[n-k])* p^(k+1)
}
return(P[N])
}
seq_dice()

seq_diceJ <- function(N=100,k=5,p=1/2){ # Just k sequence
seq_dice(N,k,p)-seq_dice(N,k+1,p)
}
seq_diceJ()
seq_diceJ(100,5,1/2)
seq_diceJ(1000000,18,1/2)
#
vsdJ=Vectorize(seq_diceJ)
N=1e5
kk=1:30
p=0.5
y=vsdJ(N,kk,0.5)
which.max(y) # 1e2:5 1e3:9 1e4:12 1e5:15 1e6:18 1e7:22 1e8:25
plot(kk,y,bty='l',pch=19,xlab='sequence',ylab='probability')
cbind(kk,y)
max(y)

Nd2mps <- function(N){ # N dice to most probable sequence
vsdJ=Vectorize(seq_diceJ)
k=min(N/2,30)
kk=1:min(k)
y=vsdJ(N,kk,0.5)
# c(seq=which.max(y),prob=max(y))
return(which.max(y))
}
Nd2mps(777)
vN2s=Vectorize(Nd2mps)
nn=seq(200,600,by=1)
y=vN2s(nn)
range(nn[y==7])
plot(nn,y,bty='l')

219132人目の素数さん2018/11/20(火) 12:14:23.98ID:cFR1wwH3
abura <- function(Amax=10,Bmax=7,Cmax=3,a=10,b=0,c=0,A=5,B=5,C=0){
# http://tpcg.io/rvz4qS
# Amax=10 ; Bmax=7 ; Cmax=3 # capacity
# t0 = c(10,0,0) # initial state
# goal = c(5,5,0) # goal
t0=c(a,b,c)
goal=c(A,B,C)

.stack <<- t(as.matrix(t0)) # stack (converted to one row matrix)
.checked <<- .stack # checked list

# pouring methods
a2b <- function(abc){ # pour from A to B
a=abc[1];b=abc[2];c=abc[3]
if(a+b > Bmax) c(a+b-Bmax,Bmax,c)
else c(0, a+b, c)
}

a2c <- function(abc){
a=abc[1];b=abc[2];c=abc[3]
if(a+c > Cmax) c(a+c-Cmax, b, Cmax)
else c(0, b, a+c)
}

b2c <- function(abc){
a=abc[1];b=abc[2];c=abc[3]
if(b+c > Cmax) c(a, b+c-Cmax,Cmax)
else c(a, 0, b+c)
}

b2a <- function(abc){
a=abc[1];b=abc[2];c=abc[3]
if(b+a > Amax) c(Amax, b+a-Amax,c)
else c(b+a, 0, c)
}

c2a <- function(abc){
a=abc[1];b=abc[2];c=abc[3]
if(c+a > Amax) c(Amax, b, c+a-Amax)
else c(c+a, b, 0)
}

c2b <- function(abc){
a=abc[1];b=abc[2];c=abc[3]
if(c+b > Bmax) c(a, Bmax, c+b-Bmax)
else c(a, c+b, 0)
}

actions = c(a2b,a2c,b2c,c2a,c2a,c2b) # all pouring method

220132人目の素数さん2018/11/20(火) 12:14:40.41ID:cFR1wwH3
is.rim <- function(row,matrix){ # is row in matrix?
any(apply(matrix,1,function(x) all(x==row))) # comparble to %in%
}

is.goal <- function(x){ # goal reached?
all(x==goal)
}

pop <- function(){ # pop LIFO
if(is.null(.stack)) return()
LIFO=.stack[1,]
if(nrow(.stack)==1) .stack <<- NULL
else .stack <<- .stack[-1,] # changed GLOBALLY
return(LIFO)
}

push <- function(rows){ # push rows at head of stack
if(is.null(rows)) invisible(NULL) # no NULL show
else .stack <<- rbind(rows , .stack) # changed GLOBELY
}

transfer <- function(x){ # return unchekcked transferred state
re=NULL
for(fun in actions){ # try all methods
v=fun(x) # drop checked state and itself
if(!is.rim(v,.checked) & !all(v==x)) re=rbind(re,fun(x))
}
uni.re=unique(re) # delete duplicated
.checked <<- rbind(uni.re,.checked) # add to .checked GLOBELY
return(uni.re)
}

state=NULL
while(!is.goal(.stack[1,])){
push(transfer(pop()))
state=rbind(state,.stack[1,])
}
return(state)
}

新着レスの表示
レスを投稿する