臨床統計もおもしろいですよ、その１ [無断転載禁止]©2ch.net
■ このスレッドは過去ログ倉庫に格納されています
0001卵の名無しさん
垢版 |
2017/05/03(水) 20:04:54.62ID:0YB5L7xG
　
　内科認定医受験の最低限の知識、
　製薬会社の示してくる臨床データ、
　論文の考察、
　論文を書くときの正当性、
　というのが、臨床統計の今までの目的の大きい部分でしたが、
　
　AI=機械学習の基本も、結局は統計学と確率に支配されます。
　そういう雑多な話をするスレです。
　
0533卵の名無しさん
垢版 |
2018/06/13(水) 08:53:38.73ID:rcG4xYvl
library(rjags)
N=50
z=40
FP=0.01
shape1=1
shape2=1
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("TP","x","theta"), n.iter=100000 )
js=as.matrix(codaSamples)
head(js)
BEST::plotPost(js[,'TP'],xlab='sensitivity')
BEST::plotPost(js[,'x'],xlab='prevalence')
BEST::plotPost(js[,'theta'],xlab='positive result',showMode = TRUE)
0534卵の名無しさん
垢版 |
2018/06/16(土) 16:23:17.07ID:V7A3qxKV
経理課の須藤は着服をやめろ！

勤務実態もないのに、グループ病院内から管理手当て（10万円）をもらうな！！！

意図的な給与操作、どうにかしろ！
0535卵の名無しさん
垢版 |
2018/06/17(日) 05:51:35.25ID:ifh2AARM
seqN <- function(N=100,K=5){
a=numeric(N)
for(i in 1:K) a[i]=2^(i-1)
for(i in K:(N-1)){
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=a[i+1]+a[i-j] # recursion formula
}
}

P0=numeric(N)
for(i in 1:N) P0[i]=a[i]/2^i # P0(n)=a(n)/2^n
P0

MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,1]=P0
head(MP);tail(MP)
MP[1,2]=1/2
for(i in (K-2):K) MP[1,i]=0

for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=1/2*MP[i,k-1]
} # Pk(n+1)=1/2*P(k-1)(n)
ret=1-apply(MP,1,sum)

ret[N]
}
seqN(100,5)
seqN(1000,10)
0536卵の名無しさん
垢版 |
2018/06/17(日) 07:58:05.56ID:ifh2AARM
## 表の出る確率がｐであるとき、N回コインを投げて　K回以上表が連続する確率
seqNp <- function(N=100,K=5,p=0.5){
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q

for(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}

P0=numeric(N)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
P0

MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
head(MP);tail(MP)
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0

for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)

ret[N]
}
0537卵の名無しさん
垢版 |
2018/06/17(日) 08:33:04.03ID:LIEnVEKd
p:表
q=1-p
# Pk(n) (k=0,1,2,3,4)を途中、5連続して表が出ていなくて
# 最後のk回は連続して表が出ている確率とする。
#
P0(1)=q
P1(1)=p
P2(1)=P3(1)=P4(1)=0
P(k+1)(n+1)=p*Pk(n)
P0(n+1)=q*{P0(n)+P1(n)+P2(n)+P3(n)+P4(n)}
=q*{P0(n)+p*P0(n-1)+p^2*P0(n-2)+p^3*P0(n-3)+p^4*P0(n-4)}

P0(n)=a(n)*p^n
# a(n+1)p^(n+1)=q*p^n{a(n)+a(n-1)+a(n-2)+a(n-3)+a(n-4)}
# a(n+1)=q/p1*(a(n)+a(n-1)+a(n-2)+a(n-3)+a(n-4))

a(n)=P0(n)/p^n
0538卵の名無しさん
垢版 |
2018/06/17(日) 09:43:49.13ID:ifh2AARM
>>532

統計くらいできるのが国立卒の普通の臨床医。

おい、ド底辺

統計処理からはおまえは

都外のド底辺シリツ医大卒と推測されたが、あってるか？
0539卵の名無しさん
垢版 |
2018/06/17(日) 09:45:16.20ID:ifh2AARM
## 表の出る確率がｐであるとき、N回コインを投げて　K回以上表が連続する確率に一般化してみた。
seqNp <- function(N=100,K=5,p=0.5){
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q
for(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}

P0=numeric(N)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
head(MP);tail(MP)
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0

for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)

ret[N]
}
0540卵の名無しさん
垢版 |
2018/06/17(日) 22:30:10.45ID:ifh2AARM
# pdfからcdfの逆関数を作ってHDIを表示させて逆関数を返す
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE){
nxx=1001
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
print(c(hdi[1],hdi[2]),digits=5)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF')
par(mfrow=c(1,1))
}
invisible(ICDF)
}
0541卵の名無しさん
垢版 |
2018/06/18(月) 20:41:20.64ID:G5tabFnb
library(rjags)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
N=100　; FP=0.01 ; shape1=1 ; shape2=1

guess.TP <- function(z){
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples(jagsModel,variable=c("TP","x"), n.iter=10000)
js=as.matrix(codaSamples)
m.TP=mean(js[,'TP'])
ci.TP=HPDinterval(as.mcmc(js[,'TP']))
m.x=mean(js[,'x'])
ci.x=HPDinterval(as.mcmc(js[,'x']))
c(m.TP=m.TP,ci.TP=ci.TP,m.x=m.x,ci.x=ci.x)
}
zz=1:20*5
re=sapply(zz,guess.TP)
head(re[,1:4])
re=as.matrix(re)
plot(zz,re['m.TP',],bty='l',ylim=c(0,1),type='n',las=1,
xlab='n : positives out of 100',ylab='sensitivity')
segments(zz,re[2,],zz,re[3,],col=8,lwd=3)
points(zz,re['m.TP',],pch=16)
0542卵の名無しさん
垢版 |
2018/06/18(月) 22:45:56.48ID:G5tabFnb
guess.TP2 <- function(z,FP){
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples(jagsModel,variable=c("TP"), n.iter=10000,thin=5)
js=as.matrix(codaSamples)
mean(js[,'TP'])
}
vG=Vectorize(guess.TP2)
n=1:20*5
FP=seq(0,0.25,length=20)
TP=outer(n,FP,vG) # wait several minutes
contour(n,FP,TP, col='navy',
xlab='n : positives out of 100',ylab='FP : 1-specificity',bty='l',nlevels=64)
points(50,0.10,pch='+',col='red',cex=1.5)
0543卵の名無しさん
垢版 |
2018/06/19(火) 22:41:00.16ID:ant1u+bV
n=5 # prime number
nn=1:(n-1)
tasu <- function(x,y) (x+y)%%n
hiku <- function(x,y) (x-y)%%n # row - col
kake <- function(x,y) (x*y)%%n
g=function(x) nn[which(x==1)]
.M=outer(nn,nn,kake)
G=apply(.M,2,g)
gyaku <- function(x) nn[which(G==(x%%n))]
waru <- function(x,y) (x*gyaku(y))%%n # row / col
waru(3,2)

xx=yy=c(0,nn)
names(xx)=paste0('x',c(0,nn))
names(yy)=paste0('y',c(0,nn))
outer(xx,yy,tasu) # x + y
outer(xx,yy,hiku) # x - y
outer(xx,yy,kake) # x * y

X=Y=nn
outer(X,Y,waru) # WRONG!!
outer(X,Y,Vectorize(waru))

a=expand.grid(X,Y)
b=matrix(mapply(waru,a[,1],a[,2]),ncol=length(X))
rownames(b)=paste0('x',nn)
colnames(b)=paste0('y',nn)
b # x / y
0544卵の名無しさん
垢版 |
2018/06/20(水) 19:25:59.38ID:myhyWcyK
rule3 <- function(n,confidence.level=0.95){
p=1/n
q=1-p # q^n.sample > 1-confidence.level
n.sample = log(1-confidence.level)/log(q)
return(n.sample)
}
n.sample=log(0.05)/log(0.999)
0545卵の名無しさん
垢版 |
2018/06/22(金) 14:47:07.52ID:ETsHYxXe
shuffle <- function(Cards){
n=length(Cards)
n1=ceiling(n/2)
n2=n-n1
C1=Cards[1:n1]
C2=Cards[(n1+1):n]
ret=NULL
for(i in 1:n1){
ret=c(ret,C1[i],C2[i])
}
ret[is.na(ret)==F]
}
x=as.character(c('A',2:10,'J','Q','K'))
cat(x,'\n') ; cat(shuffle(x))
Shuffles <- function(x){
tmp=shuffle(x)
i=1
while(!identical(x,tmp)){
tmp=shuffle(tmp)
i=i+1
}
return(i)
}

f =function(x)Shuffles(1:x)
nn=1:53
y=sapply(nn,f)
plot(nn,y,pch=16,bty='l',xlab='cards',ylab='shuffles')
cbind(nn,y)
0546卵の名無しさん
垢版 |
2018/06/22(金) 20:07:08.88ID:ETsHYxXe
http://000013.blogspot.com/2010/12/99.html

inversion <- function(x){ # 転倒数
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret)
}
x=c(4, 3, 5, 2, 1)
inversion(x)
is.even= function(x) !inversion(x)%%2
is.even(x)

prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx]) ; is.even(X)
Y=numeric(n-1)
for (i in 1:(n-1)){
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # 囚人iが見えない番号
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2]) # 偶順列になるように選択
}
all(X==Y)
}
mean(replicate(1e3,prisoner99()))
0547卵の名無しさん
垢版 |
2018/06/26(火) 07:03:32.73ID:kT/81/Gu
# http://000013.blogspot.com/2010/12/99.html
inversion <- function(x){
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret) # inversion number
}
is.even= function(x) !inversion(x)%%2 # is inverion number even?

prisoner99 <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx])
Y=numeric(n-1)
for (i in 1:(n-1)){ # select as even permutation
x1=X[-i]
x2=(1:n)[!(1:n) %in% x1] # two numbers unseen for i-th prisoner
tmp=X
tmp[i]=x2[1] ; tmp[n]=x2[2]
Y[i]=ifelse(is.even(tmp), x2[1],x2[2])
}
all(X==Y)
}
mean(replicate(1e3,prisoner99()))
0548卵の名無しさん
垢版 |
2018/06/26(火) 10:42:11.82ID:kT/81/Gu
inversion <- function(x){ #inversion number
n=length(x)
ret=numeric(n)
for(i in 1:(n-1)){
ret[i] = sum(x[i] > x[(i+1):n])
}
sum(ret)
}
is.even <- function(x) !inversion(x)%%2 # even inversion?

even.perm <- function(n=100){
indx=sample(1:n,1) # defective number
X=sample((1:n)[-indx]) # row of 99 prisoner numbers
is.even(X)
}
mean(replicate(1e3,even.perm())) # probability of even permutation
0549卵の名無しさん
垢版 |
2018/06/27(水) 14:01:31.44ID:gge/PUDl
#ある大学の学生数は５００以上１０００人以下であることはわかっている。
#無作為に２人を抽出して調べたところ
#二人とも女子学生である確率は１／２であった。
#この大学の女子学生数と男子学生数は何人か？
girlsboys <- function(g,b) g*(g-1)/(g+b)/(g+b-1)==1/2
gr=expand.grid(1:1000,1:1000)
(re=gr[which(mapply(girlsboys,gr[,1],gr[,2])),])
girlsboys(re[nrow(re),1],re[nrow(re),2])
0550卵の名無しさん
垢版 |
2018/06/27(水) 14:03:56.49ID:gge/PUDl
#　ある大学の学生数は５００以上１０００人以下であることはわかっている。
#　無作為に２人を抽出して調べたところ
#　二人とも女子学生である確率は１／２であった。
#　この大学の女子学生数と男子学生数は何人か？

girlsboys <- function(g,b) g*(g-1)/(g+b)/(g+b-1)==1/2
gr=expand.grid(1:1000,1:1000)
(re=gr[which(mapply(girlsboys,gr[,1],gr[,2])),])

#　検証
Vectorize(girlsboys)(re[,1],re[,2])
0551卵の名無しさん
垢版 |
2018/06/27(水) 17:53:09.08ID:gge/PUDl
N=120
r=10
D=c(rep(1,r),rep(0,N-r))
hiseiki <- function(m){
found=0
for(i in 1:(N-r+m)){
found=found+sample(D,1)
if(found==m) break
}
return(i)
}

re=replicate(1e4,hiseiki(3))
mean(re)
sd(re)
BEST::plotPost(re,breaks=30)
0552卵の名無しさん
垢版 |
2018/06/29(金) 07:33:02.38ID:zpGshx3p
cereals <- function(n=5){
coupons=NULL
while(!all((1:n) %in% coupons)){
coupons=append(sample(1:n,1),coupons)
}
return(length(coupons))
}

re=replicate(100,mean(replicate(1e3,cereals(5))))

> mean(re)
[1] 11.43503
0553卵の名無しさん
垢版 |
2018/06/29(金) 08:13:57.68ID:zpGshx3p
p=4:1
cereals <- function(p){
n=length(p)
coupons=NULL
while(!all((1:n) %in% coupons)){
coupons=append(sample(1:n,1,p=p),coupons)
}
return(length(coupons))
}

mean(replicate(1e3,cereals(p)))

re=replicate(100,mean(replicate(1e3,cereals(p))))

mean(re)
0554卵の名無しさん
垢版 |
2018/06/29(金) 16:00:32.65ID:zpGshx3p
blood.type <- function(p,need){
n=length(p)
ABO=NULL
enough <- function(x){
pool=numeric(n)
for(i in 1:n) pool[i]=sum(ABO==i)
all(pool >= need)
}
while(!enough(ABO)){
ABO=append(sample(1:n,1,p=p),ABO)
}
return(length(ABO))
}
p=4:1
need=c(10,10,5,2)
re=replicate(1e4,blood.type(p,need))
BEST::plotPost(re)
0555卵の名無しさん
垢版 |
2018/06/30(土) 05:47:26.70ID:BDmYWstD
BT <- function (a,b,c,d) 1/a + 1/b + 1/c + 1/d - 1/(a+b) - 1/(a+c) - 1/(b+c) - 1/(a+d) - 1/(b+d) - 1/(c+d) + 1/(a+b+c) + 1/(d+a+b) + 1/(c+d+a) + 1/(b+c+d) - 1/(a+b+c+d)
a=4
b=3
c=2
d=1
s =a+b+c+d
a=a/s
b=b/s
c=c/s
d=d/s
BT(a,b,c,d)
0556卵の名無しさん
垢版 |
2018/06/30(土) 13:48:31.26ID:j2CU1Lw0
even.tally <- function(a=3 , b=2){
idx=combn(1:(a+b),a)
n=ncol(idx)
mat=matrix(0,nrow=n,ncol=a+b)
for(i in 1:n) mat[i,idx[,i]]=1
tally <- function (x) any(cumsum(x)==cumsum(1-x))
mean(apply (mat,1,tally))
}
even.tally()
even.tally(5,10)
0557卵の名無しさん
垢版 |
2018/06/30(土) 15:09:55.32ID:j2CU1Lw0
a=750 ; b=250
v=c(rep(1,a),rep(0,b))

f <- function(v){
x=sample(v)
any(cumsum(x)==cumsum(1-x))
}

mean(replicate(1e5,f(v)))
0558卵の名無しさん
垢版 |
2018/07/01(日) 10:28:46.39ID:F3KIhVUK
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x)
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}

#
birth(100, 3)

vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
0559卵の名無しさん
垢版 |
2018/07/01(日) 11:53:44.57ID:F3KIhVUK
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x))
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}

#
birth(100, 3)

vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
0560卵の名無しさん
垢版 |
2018/07/01(日) 12:03:30.64ID:F3KIhVUK
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x)
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}

#
birth(100, 3)

vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
0561卵の名無しさん
垢版 |
2018/07/01(日) 12:08:34.56ID:F3KIhVUK
インフルエンザの迅速キットは特異度は高いが感度は検査時期によって左右される。
ある診断キットが開発されたとする。
このキットは特異度は９９％と良好であったが、
感度については確かな情報がない。
事前確率分布として一様分布を仮定する。
５０人を無作為抽出してこの診断キットで診断したところ４０人が陽性であった。
この母集団の有病率の期待値と９５％信用区間はいくらか？
またこの診断キットの感度の期待値と９５％信用区間はいくらか



暇つぶしにこれをMCMCを使わずに解く方法を考えていた。
偽陽性率FP=0.01として
陽性確率p=TP*x+(1-x)*FP
尤度が50C40*p^40*(1-p)^10
TPは一様分布なので積分消去して
確率密度関数に比例する関数を作ってarea under the curveで割って確率密度関数化したのち積分して累積密度関数をつくる。この累積密度関数の逆関数を作って95%区間が最短になる区間を計算すれば信頼区間が算出できる。


この結果がstanでのシミュレーションの結果と一致すればよし。
0562卵の名無しさん
垢版 |
2018/07/01(日) 22:00:46.73ID:lSk1MYzX
# # choose(n,r) == gamma(n+1) / (gamma(r+1) * gamma(n-r+1))
same.birthday <- function(n) 1-choose(365+97/400,n)*factorial(n)/(365+97/400)^n
plot(x,y,bty='l',xlab='subjects',ylab='probability')
curve(same.birthday(x),add = TRUE)
abline(h=0.5,col=8)
same.birthday(22:23)
0563卵の名無しさん
垢版 |
2018/07/02(月) 01:26:48.74ID:pp47QgIN
library(rjags)
N=50
z=40
FP=0.01
shape1=1
shape2=1
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("TP","x","theta"), n.iter=100000 )
js=as.matrix(codaSamples)
head(js)
BEST::plotPost(js[,'TP'],xlab='sensitivity')
BEST::plotPost(js[,'x'],xlab='prevalence')
BEST::plotPost(js[,'theta'],xlab='positive result',showMode = TRUE)
0564卵の名無しさん
垢版 |
2018/07/02(月) 01:28:46.92ID:pp47QgIN
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
real<lower=0,upper=1> x; //prevalence
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=TP*x+FP*(1-x);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2); // T[0.5,];
x ~ beta(shape1,shape2);
}

')
0565卵の名無しさん
垢版 |
2018/07/02(月) 01:29:10.57ID:pp47QgIN
# model=stan_model(model_code = stanString)
# saveRDS(model,'quick_kit.rds')
model=readRDS('quick_kit.rds')
fit=sampling(model,data=data,iter=10000)
print(fit,digits=3,probs=c(.025,.50,.975))
stan_trace(fit)
# stan_diag(fit)
stan_ac(fit)
stan_dens(fit,separate_chains = TRUE)
stan_hist(fit,fill='skyblue',bins=15,pars=c('x','TP'))
ms=rstan::extract(fit)
BEST::plotPost(ms$TP,showMode = TRUE,xlab='TP')
BEST::plotPost(ms$x,showMode = FALSE,xlab='prevalence',col=sample(colours(),1))
0566卵の名無しさん
垢版 |
2018/07/02(月) 21:40:12.88ID:tMXA02ZR
escape.cliff <- function(p=2/3,k=1000){
pos=1
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e2,escape.cliff()))
0567卵の名無しさん
垢版 |
2018/07/02(月) 22:58:14.27ID:tMXA02ZR
escape.cliff <- function(pos=1,p=2/3,k=1000){
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e3,escape.cliff(10,0.5,100)))

#
totter <- function(p=2/3,pos=1,pos0=NULL,k=1e3){
if(is.null(pos0)) pos0=pos
for(i in 1:k) {
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(pos>pos0)
}
mean(replicate(1e3,totter(0.5,10,5 )))
0568卵の名無しさん
垢版 |
2018/07/03(火) 07:46:06.00ID:22JZXDLY
exam <- function(p=0.5,hit=0,money=5,k=5){
while(hit<k){
money=money-1
if(money==0) return(FALSE)
shoot=sample(c(1,0),1,prob=c(p,1-p))
if(shoot){
money=money+1
hit=hit +1
}
}
return(TRUE)
}
mean(replicate(1e5, exam()))
0569卵の名無しさん
垢版 |
2018/07/03(火) 08:47:34.93ID:ltuSOSv2
ryunen <- function(p=0.6 ,money=10,grade=0,ryu=0){
while(grade<6){
test=rbinom(1,1,p)
if(test) grade=grade+1
else ryu=ryu+1
if(ryu==2){
money=money-5
ryu=ryu-1
}
if(money<=0) return(FALSE)
}
return(grade==6)
}
mean(replicate(1e5,ryunen()))
0570卵の名無しさん
垢版 |
2018/07/03(火) 11:46:57.62ID:22JZXDLY
escape.cliff <- function(pos=1,p=2/3,k=1000){
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e5,escape.cliff(1,2/3,3)))
4/7
0571卵の名無しさん
垢版 |
2018/07/03(火) 13:41:51.20ID:22JZXDLY
gmbl <- function(money=20,p=18/38){
while(0<money & money <40){
money=money+sample(c(1,-1),1,p=c(p,1-p))
}
return(money==40)
}
mean(replicate(1e3, gmbl()))
0572卵の名無しさん
垢版 |
2018/07/03(火) 14:41:00.85ID:22JZXDLY
x=c(rep(1,4),rep(0,48))
mean(replicate(1e5,which(sample(x)==1)[1]))


y=numeric(48)
for(i in 1:48) y[i]=i*choose(48,i-1)*factorial(i-1)*4*factorial(52-i)/factorial(52)
sum(y)
0573卵の名無しさん
垢版 |
2018/07/03(火) 14:42:24.05ID:22JZXDLY
(a) A railroad numbers its locomotives in order, 1, 2, . . . , N. One day you see a locomotive and its number is 60.
Guess how many locomotives the company has.
(b) You have looked at 5 locomotives and the largest number observed is 60.
Again guess how many locomotives the company has.
0574卵の名無しさん
垢版 |
2018/07/03(火) 16:06:53.54ID:22JZXDLY
n=5
m=60
N=n:100
loco <- function(x){
max(sample(1:x,n))
}
vloco=Vectorize(loco)
loco.sim <- function() N[which(vloco(N)==m)]
locomotives=unlist(replicate(1e4,loco.sim()))
BEST::plotPost(locomotives,breaks=30)
HDInterval::hdi(locomotives)
0575卵の名無しさん
垢版 |
2018/07/03(火) 18:26:54.88ID:22JZXDLY
n=5
library(gtools)
y=permutations(n,n,1:n)
f=function(x) sum(x==1:n)
sum(apply(y,1,f))/length(y)
1/n
0576卵の名無しさん
垢版 |
2018/07/03(火) 18:39:40.14ID:22JZXDLY
f = function(n){
sum(sample(1:n)==1:n)/n
}
g=function(n)mean(replicate(1e4,f(n)))
vg=Vectorize(g)
x=1:100
plot(x,vg(x))
curve(1/x,add=T)
0577卵の名無しさん
垢版 |
2018/07/03(火) 21:26:39.83ID:22JZXDLY
gr=expand.grid(1:9,1:9,1:9)
u=10^(2:0)
f=function(a,b,c,x=1776) sum(c(a,b,c)*u+c(b,c,a)*u+c(c,a,b)*u)-x
# when A>B>C
idx=which(mapply(f,gr[,1],gr[,2],gr[,3])==0 & gr[,1]>gr[,2] & gr[,2]>gr[,3])
gr[idx,]
# otherwise
idx2=which(mapply(f,gr[,1],gr[,2],gr[,3])==0)
length (idx2)
gr[idx2,]
0578卵の名無しさん
垢版 |
2018/07/04(水) 00:37:12.00ID:sDwIwpA3
n=5
m=60
N=n:100
loco <- function(x){
max(sample(1:x,n))
}
vloco=Vectorize(loco)
loco.sim <- function() N[which(vloco(N)==m)]
locomotives=unlist(replicate(1e4,loco.sim()))
BEST::plotPost(locomotives,breaks=30)

HDInterval::hdi(locomotives)
n=60:100
pmf=choose(59,4)/choose(n,5) #Pr(max=60|n)
pdf=pmf/sum (pmf)
sum( n*pdf) #E(n)
lines(n,pdf)
0579卵の名無しさん
垢版 |
2018/07/05(木) 10:51:42.59ID:9YHvFei/
own <- function(n){
mean(replicate(1e5,sum(sample(1:n)-1:n==0)))
}

own(100)
0580卵の名無しさん
垢版 |
2018/07/05(木) 10:54:03.43ID:9YHvFei/
r1 <- function(x){ # rotate by one bead
n=length(x)
y=numeric(n)
y[1]=x[n]
for(i in 1:(n-1)){
y[i+1]=x[i]
}
return(y)
}

rn <- function(x){ # every rotation
n=length(x)
mat=matrix(rep(NA,n^2),ncol=n)
mat[,1]=x
for(i in 1:(n-1)){
mat[,i+1]=r1(mat[,i])
}
return(t(mat))
}

same <- function(x,y){
if(sum(x)!=sum(y)) return(FALSE)
f=function(a,b=y){ # is equal to y
all(a==b)
}
mat=rbind(rn(x),rn(rev(x))) # with symmetric conversion
any(apply(mat,1,f))
}
0581卵の名無しさん
垢版 |
2018/07/05(木) 10:55:44.62ID:9YHvFei/
dec2bin <- function(num, digit=0){ # decimal to 0,1 vector
if(num <= 0 && digit <= 0){
return(NULL)
}else{
return(append(Recall(num%/%2,digit-1), num%%2))
}
}
vd2b=Vectorize(dec2bin)
#
bracelett <- function(n){
mat=t(vd2b(0:(2^n-1),n)) # make all permutation of beads
# head(mat)
ret=list() # list of the same bracelett
for(i in 1:2^n){
ret[[i]]=which(apply(mat,1,function(z)same(z,mat[i,])))
}
# head(ret) ; table(unlist(ret))
del=NULL
for(i in 1:2^n) del=append(del,ret[[i]][-1])
2^n - length(unique(del))
}
0582卵の名無しさん
垢版 |
2018/07/05(木) 14:59:01.13ID:bn3tqJqU
> library(gtools)
> swap <- function(n){
+ perm=permutations(n,n,1:n)
+ mean(apply(perm,1,function(x)sum(x==1:n)))
+ }
> swap(9)
[1] 1
> swap.sim <- function(n,k=1e5){
+ mean(replicate(k,sum(sample(1:n)==1:n)))
+ }
> swap.sim(9)
[1] 1.0004
0583卵の名無しさん
垢版 |
2018/07/05(木) 20:18:51.37ID:bn3tqJqU
swap2 <- function(n){ # k*nCk*(1/n)^k*(1-1/n)^(n-k) k=0:n
re=numeric(n)
for(i in 0:n){
p=1/n
re[i]=i*choose(n,i)*p^i*(1-p)^(n-i)
}
sum(re)
}
swap2(7)
#
swap3 <- function(n){
re=numeric(n)
for(i in 0:n){
p=1/n
re[i]=i*dbinom(i,n,p)
}
sum(re)
}
swap3(7)
0584卵の名無しさん
垢版 |
2018/07/06(金) 09:11:30.20ID:UEcW6fma
f012 <- function(n){
args=list()
length(args)=n-1
args[[1]]=args[[n-1]]=1:2
for(i in 2:(n-2)){
args[[i]]=0:2
}
gr=do.call(expand.grid,args)
gr=as.matrix(gr)
ret=numeric()
for(i in 1:nrow(gr)){
ret[i]=all(diff(gr[i,])!=0)
}
sum(ret)
}
0585卵の名無しさん
垢版 |
2018/07/21(土) 16:31:14.15ID:V1Wm3iKf
/* SEND+MORE=MONEYの覆面暗算をC言語で解いてみる。*/
#include<stdio.h> /* SEND+MORE=MONEY */
int compare_int(const void *a, const void *b)
{return *(int*)a - *(int*)b;}
int i,j;
int unique(int num[8]){
qsort(num,8,sizeof(int),compare_int);
for(i=0;i<8;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;
}
main(){
int S,E,N,D,M,O,R,Y;
for(S = 1; S < 10; S++){
for(E = 0; E < 10; E++){
for(N = 0; N < 10; N++){
for(D = 0; D < 10; D++){
for(M = 1; M < 10; M++){
for(O = 0; O < 10; O++){
for(R = 0; R < 10; R++){
for(Y = 0; Y < 10; Y++){
if(S*1000+E*100+N*10+D+M*1000+O*100+R*10+E==M*10000+O*1000+N*100+E*10+Y){
int num[]={S,E,N,D,M,O,R,Y};
if(unique(num)==1){
printf("%d%d%d%d+%d%d%d%d=%d%d%d%d%d\n", S,E,N,D,M,O,R,E,M,O,N,E,Y);
}}}}}}}}}}}
C:\MinGW>gcc sendmore.c -o money
C:\MinGW>money
9567+1085=10652
0586卵の名無しさん
垢版 |
2018/07/22(日) 11:46:07.83ID:tW9s1ZUi
#include<stdio.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;}
int unique(int num[]){
int i,j,n=10;
qsort(num,n,sizeof(int),compare_int);
for(i=0;i<n;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;}
main(){ /* ド底辺＋私立医＝裏口馬鹿 */
int n=1,A,B,C,D,E,F,G,H,I,J;
for(A = 1; A < 10; A++){
for(B = 0; B < 10; B++){
for(C = 0; C < 10; C++){
for(D = 1; D < 10; D++){
for(E = 0; E < 10; E++){
for(F = 0; F < 10; F++){
for(G = 1; G < 10; G++){
for(H = 0; H < 10; H++){
for(I = 0; I < 10; I++){
for(J = 0; J < 10; J++){
if(A*100+B*10+C +D*100+E*10+F==G*1000+H*100+I*10+J){
int num[]={A,B,C,D,E,F,G,H,I,J};
if(unique(num)==1){
printf("%2d: %d%d%d + %d%d%d = %d%d%d%d\n", n,A,B,C,D,E,F,G,H,I,J);
n++;
}}}}} }}} }}} }}
0587卵の名無しさん
垢版 |
2018/07/23(月) 16:11:37.39ID:B1Q5jsGg
# ２つの整数があります。
# それらをたしてできた数は、
# 十の位と一の位の数字が等しい２けたの整数になり、
# それらをかけてできた数は、
# 百の位、十の位、一の位が等しい３けたの整数になりました。
# このような２つの整数の組をあるだけ答えなさい。
gr=expand.grid(1:99,1:99)
f <- function(x,y){ all(x>=y,c((x+y)%%10==(x+y)%/%10,(x*y)%/%100==((x*y)%/%10)%%10,(x*y)%/%100==(x*y)%%10))}
gr[which(mapply(f,gr[,1],gr[,2])==TRUE),]
0588卵の名無しさん
垢版 |
2018/07/23(月) 17:34:03.68ID:B1Q5jsGg
# 実行速度遅すぎて実用性なし

# Consider a function which, for a given whole number n,
# returns the number of ones required when writing out all numbers between 0 and n.
# For example, f(13)=6. Notice that f(1)=1.
# What is the next largest number n such that f(n)=n.

ｆ０　＜−　ｆｕｎｃｔｉｏｎ（ｎ）｛
　　ｙ＝ａｓ．ｃｈａｒａｃｔｅｒ（ｎ）
　　ｎｃ＝ｎｃｈａｒ（ｙ）
　　ｚ＝ＮＵＬＬ
　　ｆｏｒ（ｉ　ｉｎ　１：ｎｃ）　ｚ［ｉ］　＜−　ｓｕｂｓｔｒ（ｙ，ｉ，ｉ）＝＝１
　　ｒｅｔｕｒｎ（ｓｕｍ（ｚ））
｝
ｆ　＜−　ｆｕｎｃｔｉｏｎ（ｎ）｛
　　ｒｅ＝ｎｕｍｅｒｉｃ（）
　　ｒｅ［１］＝１
　　ｆｏｒ（ｉ　ｉｎ　２：ｎ）｛
　　　　ｒｅ［ｉ］＝ｒｅ［ｉ−１］＋ｆ０（ｉ）
　　｝
　　ｒｅｔｕｒｎ（ｒｅ［ｎ］）
｝
ｇ　＜−　ｆｕｎｃｔｉｏｎ（ｎ）　　ｎ　−　ｆ（ｎ）
ｉ＝１９９９７９
ｗｈｉｌｅ（ｇ（ｉ）！＝０）｛
　　ｉ＜−ｉ＋１
｝
ｉ
0589卵の名無しさん
垢版 |
2018/07/23(月) 22:11:59.00ID:B1Q5jsGg
6!=(5+1)5!=5*5!+5!
=5*5!+(4+1)*4!
=5*5!+4*4!+4!
=5*5!+4*4!+(3+1)*3!
=5*5!+4*4!+3*3!+3!
=5*5!+4*4!+3*3!+(2+1)*2!
=5*5!+4*4!+3*3!+2*2!+2!
=5*5!+4*4!+3*3!+2*2!+1*1!+1

factorial(n+1)= sum(factoral(n:1)*(n:1))+1
0590卵の名無しさん
垢版 |
2018/07/24(火) 09:34:20.21ID:RZUDxWdB
Ｂ，Ｃ，Ｄ，Ｅ，Ｆが０〜９の数字（同じ数字であってもよい）で
６！＊Ｂ＋５！＊Ｃ＋４！＊Ｄ＋３！＊Ｅ＋２！＊Ｆ＋１！＊Ｇ＝５５５５
が成立するときＢ＋Ｃ＋Ｄ＋Ｅ＋Ｆ＋Ｇの最小となるＢ〜Ｇの組み合わせを求めよ。

n=5555
f <- function(B,C,D,E,F,G) sum(factorial(6:1)*c(B,C,D,E,F,G))-n
max=n%/%factorial(6)
gr=expand.grid(0:max,0:6,0:5,0:4,0:3,0:2)
ret=mapply(f, gr[,1],gr[,2],gr[,3],gr[,4],gr[,5],gr[,6])
bg=apply(gr[which(ret==0),],1,sum)
(min.sum=bg[which.min(bg)])
gr[names(min.sum),]

n%/%factorial(6)
n%%factorial(6)%/%factorial(5)
n%%factorial(6)%%factorial(5)%/%factorial(4)
n%%factorial(6)%%factorial(5)%%factorial(4)%/%factorial(3)
n%%factorial(6)%%factorial(5)%%factorial(4)%%factorial(3)%/%factorial(2)
n%%factorial(6)%%factorial(5)%%factorial(4)%%factorial(3)%%factorial(2)%/%factorial(1)
0591卵の名無しさん
垢版 |
2018/07/24(火) 16:06:02.90ID:Cp/BkN3V
#include<stdio.h>
#include<string.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;
}

main( int argc, char *argv[] ){
char N = '0'+ atoi(argv[1]);
int a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t;
for(a = 0; a < 10; a++){
for(b = 0; b < 10; b++){
for(c = 0; c < 10; c++){
for(d = 0; d < 10; d++){
for(e = 0; e < 10; e++){
for(f = 0; f < 10; f++){
for(g = 0; g < 10; g++){
for(h = 0; h < 10; h++){
for(i = 0; i < 10; i++){
for(j = 0; j < 10; j++){
for(k = 0; k < 10; k++){
for(l = 0; l < 10; l++){
0592卵の名無しさん
垢版 |
2018/07/24(火) 16:07:14.87ID:Cp/BkN3V
for(m = 0; m < 10; m++){
for(n = 0; n < 10; n++){
for(o = 0; o < 10; o++){
for(p = 0; p < 10; p++){
for(q = 0; q < 10; q++){
for(r = 0; r < 10; r++){
for(s = 0; s < 10; s++){
for(t = 0; t < 10; t++){
int num[]={a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t};
qsort(num,20,sizeof(int),compare_int);
char NUM[20];
int i;
for(i=0;i<20;i++){
NUM[i]='0'+ num[i];
}
char N4[4]={N,N,N,N};
char N5[5]={N,N,N,N,N};
if(strstr(NUM,N4)!=NULL & strstr(NUM,N5)==NULL &
(100*a+10*b+c)*f == 100*g+10*h+i &
(100*a+10*b+c)*e == 100*j+10*k+l &
(100*a+10*b+c)*d == 100*m+10*n+o &
(100*m+10*n+o)*100 + (100*j+10*k+l)*10 + 100*g+10*h+i == p*10000+q*1000+r*100+s*10+t
){
printf("abc def ghi jkl mno pqrst = %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d%d%d\n",
a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t);
}
}}}}} }}}}} }}}}} }}}}}
}
0593卵の名無しさん
垢版 |
2018/07/28(土) 06:54:00.42ID:Y3H/zTDn
options(scipen = 10)
fibon <- function(N){
f=numeric(N)
f[1]=1
f[2]=1
if(N>2){for(i in 3:N){
f[i]=f[i-2]+f[i-1]
}}
return(f[N])
}
Vectorize(fibon)(1:15)
fibon(50)
fibon(100)# wrong!
0595卵の名無しさん
垢版 |
2018/08/01(水) 09:03:07.74ID:PtTCNKBk
?- L= [1,2,3,4,5,6,7,8,9],c([_,_,_],X,L),c(Y,[_,_,_],X).
L = [1,2,3,4,5,6,7,8,9],
X = [4,5,6,7,8,9],
Y = [4,5,6] ;
No
0596卵の名無しさん
垢版 |
2018/08/01(水) 21:24:01.69ID:nhdpq/Mi
last([Last],Last).
last([_|Rest],Last) :- last(Rest,Last).
% last([1,2,3,4,5],Last).
0597卵の名無しさん
垢版 |
2018/08/04(土) 20:38:03.88ID:dITWU8BY
日本人の血液型はＡ,Ｏ,Ｂ,ＡＢの比率が４：３：２：１であるという。
それぞれの血液型の人を最低でも各々４、３、２、１人集めるためには必要な人数の期待値はいくらか？
0598卵の名無しさん
垢版 |
2018/08/04(土) 21:45:28.37ID:Tkj7u78K
>>597
blood.type <- function(p,need){
n=length(p)
ABO=NULL
enough <- function(x){
pool=numeric(n)
for(i in 1:n) pool[i]=sum(ABO==i)
all(pool >= need)
}
while(!enough(ABO)){
ABO=append(sample(1:n,1,p=p),ABO)
}
return(length(ABO))
}
p=4:1
need=c(4,3,2,1)
re=replicate(1e5,blood.type(p,need))
mean(re)
0599卵の名無しさん
垢版 |
2018/08/04(土) 21:52:26.83ID:Tkj7u78K
blood.samples <- function(p=c(4,3,2,1),need=c(1,1,1,1),k=1e5){
mean(replicate(k,blood.type(p,need)))
}
blood.samples()
blood.samples(need=4:1,k=1e4)
0600卵の名無しさん
垢版 |
2018/08/05(日) 07:20:33.69ID:aHWo4FJr
ｓｕｂｓｅｔ（［］，［］）．
ｓｕｂｓｅｔ（［Ｆｉｒｓｔ｜Ｒｅｓｔ］，［Ｆｉｒｓｔ｜Ｓｕｂ１］）　：−　ｓｕｂｓｅｔ（Ｒｅｓｔ，Ｓｕｂ１）．
ｓｕｂｓｅｔ（［＿｜Ｒｅｓｔ］，Ｓｕｂ２）　：−　ｓｕｂｓｅｔ（Ｒｅｓｔ，Ｓｕｂ２）．


？−　ｓｕｂｓｅｔ（［ド底辺，特殊，シリツ医大，イカサマ入試，裏口，馬鹿］，＿恥），　ｗｒｉｔｅｌｎ（＿恥），ｆａｉｌ．
0602卵の名無しさん
垢版 |
2018/08/05(日) 20:30:16.90ID:aHWo4FJr
#https://i.imgur.com/Y81vbNp.jpg
N=5
library(gtools)
perm=permutations(n=2,r=N,v=c(1,5), repeats.allowed = TRUE)
move=t(apply(perm,1,cumsum))
p0=0
P=(p0+move)%%6
q0=2
Q=(q0+move)%%6

is.fe <- function(x,y){ # is first encounter after N sec?
all(x[-N]!=y[-N]) & x[N]==y[N]
}
re=NULL
for(i in 1:nrow(P)){
for(j in 1:nrow(Q)){
re=append(re,is.fe(P[i,],Q[j,]))
}
}
(Ans=mean(re)) ; 3^(N-1)/(2^N)^2
cat(Ans*(2^N)^2,'/',(2^N)^2)
0603卵の名無しさん
垢版 |
2018/08/06(月) 07:05:55.55ID:29wsmNi+
is.pe <- function(x,y,M) { # encounter M or more than M times within N sec.
sum(x==y) >= M
}
M=2
re=NULL
for(i in 1:nrow(P)){
for(j in 1:nrow(Q)){
re=append(re,is.pe(P[i,],Q[j,],M))
}
}
(Ans=mean(re))
cat(Ans*(2^N)^2,'/',(2^N)^2)
0605卵の名無しさん
垢版 |
2018/08/08(水) 19:31:00.52ID:OkpyQ+1n
hanoi <- function(n,from='A',via='B',to='C'){
if(n >= 1){
Recall(n-1,from,to,via)
cat('move',n,'from',from,' to ',to,'\n')
Recall(n-1,via,from,to)
}
}
hanoi(4)
0606卵の名無しさん
垢版 |
2018/08/09(木) 05:35:20.13ID:mD+P4tQf
掛け算を再帰関数で定義する。

VAT <- function(n,m=1.08){
if(n==0) 0
else Recall(n-1) + m
}
VAT(250)
0607卵の名無しさん
垢版 |
2018/08/09(木) 06:51:20.87ID:mD+P4tQf
総和を再帰関数で書く

sup <- function(v){
if (length(v)==0) 0
else v[1] + Recall(v[-1])
}
sup(1:10)
0608卵の名無しさん
垢版 |
2018/08/09(木) 10:57:33.87ID:JutZ+/A9
# Ackermann
> A <- function(m,n){
+ if(m==0) return(n+1)
+ if(n==0) return(Recall(m-1,1))
+ else return(Recall(m-1,Recall(m,n-1)))
+ }
>
> A(2,4)
[1] 11
> A(3,4)
[1] 125
> A(4,1)
Error: evaluation nested too deeply: infinite recursion / options(expressions=)?
>
0609卵の名無しさん
垢版 |
2018/08/09(木) 18:49:05.66ID:JutZ+/A9
累積和

# cumsum with for-loop
cumsumL <- function(v){
n=length(v)
re=numeric(n)
re[1]=v[1]
for(i in 1:(n-1)) re[i+1]=re[i]+v[i+1]
re
}

# cumsum with recursive call
cumsumR <- function(v,res=NULL,i=1){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
}
0610卵の名無しさん
垢版 |
2018/08/09(木) 20:12:53.52ID:JutZ+/A9
# 10進法をN進法でdigit桁表示する

dec2n <- function(num, N = 2, digit = 0){ # decimal to 0,1,..,n-1 vector
if(num==0 & digit==0) return(0)
if(num <= 0 & digit <= 0) return()
else{
return(append(Recall(num%/%N, N ,digit-1), num%%N))
}
}

> dec2n(0)
[1] 0
> dec2n(11,digit=5)
[1] 0 0 1 0 1 1
> dec2n(9,N=5,digit=3)
[1] 0 0 1 4
> dec2n(1000,N=16)
[1] 3 14 8
0611卵の名無しさん
垢版 |
2018/08/09(木) 21:28:39.07ID:JutZ+/A9
>>610　(degugged)
dec2n <- function(num, N = 2, digit = 1){ # decimal to 0,1,..,n-1 vector
if(num <= 0 & digit <= 0) return()
else{
return(append(Recall(num%/%N, N ,digit-1), num%%N))
}
}
dec2n(0)
dec2n(11,digit=5)
dec2n(9,N=5,digit=3)
dec2n(1000,N=16)

#
dec2hex <- function(x){ # decimal to hexa
hex=c(0:9,letters[1:6])
n=length(x)
re=numeric(n)
for(i in 1:n){
if(x[i]==0) re[i]='0'
else re[i]=hex[x[i]+1]
}
cat(re,'\n')
}

dec2hex(c(1,0,15))
dec2hex(dec2n(1000,16))

dec2sexa <- function(x) dec2hex(dec2n(x,16))
0612卵の名無しさん
垢版 |
2018/08/14(火) 14:27:35.01ID:Z2jjlChF
draft

X <- function(n,red=10,white=90){
rw=red+white
total=martix(0,nrow=n,ncol=2^n)
p=martix(0,nrow=n,ncol=2^n)
total[1,1:2]=c(rw-1,rw+1)
p[1,1:2]=c(white/rw,red/rw)
if(n > 1){
for(i in 1:n){
li=2^i
total[i+1,1:(2*li)]=c(total[i,1:li]-1, total[i,1:li]+1)
p[i+1,1:2*li)]=c(p[i,1:li]*(total[i,1:li]-red)/total[i,1:li], p[i,1:li]*(red)/total[i,1:li])
}
}
return(sum(p[n,]*total[n,]))
}
0613卵の名無しさん
垢版 |
2018/08/15(水) 08:08:12.45ID:SmB+loM1
rm(list=ls())

X = function(n,r=10,w=90){ #n:試行数 r: 赤玉数 w:白玉数
rw=r+w # 試行前総玉数
J=rw+n # 総玉数の上限
# s[i,j] i回試行後に総数がj個である確率の行列
s=matrix(0,nrow=n,ncol=J)
s[1,rw-1]=w/rw ; s[1,rw+1]=r/rw # 1回試行後
if(n > 1){
for(i in 2:n){
for(j in r:J){ # jはr未満にはならない
# if(j==1) s[i,j] = s[i-1,j+1]*(j+1-r)/(j+1)
if(j==J) s[i,j] = s[i-1,j-1] * r/(j-1)
else s[i,j] = s[i-1,j-1] * r/(j-1) + s[i-1,j+1]*(j+1-r)/(j+1)
}
}
}
total=sum((r:J)*s[n,r:J])# n回試行後総数の期待値
white=total-r
return(c(total=total,white=white))
}


> vX=Vectorize(X)
> vX(1000:1010)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
total 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5
white 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5
0614卵の名無しさん
垢版 |
2018/08/15(水) 12:17:01.96ID:eDY+Gzxn
# シオマネキ
# http://i.imgur.com/uBTnujq.png
H=(sqrt(5+2*sqrt(5)))/2
A=0.5+0i
B=A+r

shiomaneki <- function(x,y){
xy=x+y*1i
(abs(xy-A)+abs(xy-B)+x)*2 + 1
}
x=seq(0,Re(B),length=100)
y=seq(0,H,length=100)
z=outer(x,y,shiomaneki)
contour(x,y,z,nlevels=20,lwd=2,bty='l')
x11() ;persp(x,y,z, theta=35, lty=3,col="lightblue",xlab='x',ylab='y',zlab='z',
ticktype='detailed',shade=0.75,phi=30,ltheta=-10,border=TRUE)
rgl::persp3d(x,y,z,col='blue')

sio <- function(xy) shiomaneki(xy[1],xy[2])
optim(c(0.4,0.4),sio,method='Nelder-Mead')
0615卵の名無しさん
垢版 |
2018/08/15(水) 12:17:36.49ID:eDY+Gzxn
# カブトガニ
# https://i.imgur.com/VhC8cug.png
H=(sqrt(5+2*sqrt(5)))/2
A=0.5+0i
B=A+r
C=0+H*1i
abs(C-B)
kabutogani <- function(y1,y2){
K=0+y1*1i
G=0.5+y2*1i
abs(K-C) + 2*abs(K-G) + 2*abs(G-A) +2*abs(G-B)
}
y1=seq(0,H,length=100)
y2=seq(0,H,length=100)
z=outer(y1,y2,kabutogani)
contour(y1,y2,z,nlevel=100,lwd=2,xlim=c(0.5,1.5),ylim=c(0,1.2),bty='l')

kab <- function(y1y2) kabutogani(y1y2[1],y1y2[2])
optim(c(1,0.5),kab, method = 'Nelder-Mead')
0616卵の名無しさん
垢版 |
2018/08/15(水) 15:14:03.13ID:eDY+Gzxn
DOP <- function(n,print=FALSE){ # diagonal length of regular plygon of size length 1
q=2*pi/n
r=cos(q)+1i*sin(q)
p=numeric(n+1)
for(i in 1:(n+1)) p[i]=r^(i-1)
D=NULL
if(n>3){
for(j in 3:(2+ceiling(n-3)/2)){
D= append(D,abs(p[1]-p[j]))
}
}
if(print){
plot(p,type='l',bty='l',axes=FALSE,ann=FALSE,lwd=2)
for(i in 3:(n-1)) segments(1,0,Re(p[i]),Im(p[i]),col='gray')
}
return(D)
}

DOP(17,p=T)
0617卵の名無しさん
垢版 |
2018/08/15(水) 22:04:58.96ID:SmB+loM1
DOP <- function(n,print=FALSE){ # diagonal length of regular plygon of size length 1
q=2*pi/n
r=cos(q)+1i*sin(q)
p=numeric(n+1)
for(i in 1:(n+1)) p[i]=r^(i-1)
D=NULL
if(n>3){
for(j in 3:(2+ceiling(n-3)/2)){
D= append(D,abs(p[1]-p[j])/abs(p[1]-p[2]))
}
}
if(print){
plot(p,type='l',bty='l',axes=FALSE,ann=FALSE,lwd=2)
for(i in 3:(n-1)) segments(1,0,Re(p[i]),Im(p[i]),col='gray')
}
return(D)
}

DOP(17,p=T)
0618卵の名無しさん
垢版 |
2018/08/15(水) 23:51:51.61ID:eDY+Gzxn
> (sin(2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - (cos(4*2*pi/7)*sin(2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - sin(4*2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) + (cos(2*pi/7)* sin(4*2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2))
[1] 2.190643
0619卵の名無しさん
垢版 |
2018/08/16(木) 22:05:33.14ID:8kKe3yXf
rm(list=ls())
graphics.off()
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,lwd=1,asp=1)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),1:n,cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
0620卵の名無しさん
垢版 |
2018/08/17(金) 09:11:30.66ID:qXmvzl8y
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),paste('p',1:n),cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),col=2,lwd=2,...)
pt <- function(x,y,...) text(Re(x),Im(x), y,...)

kabutogani3 <- function(xl,yl,xc,yc,xr,yr){
L=xl+yl*1i
C=xc+yc*1i
R=xr+yr*1i
abs(p[3]-C)+abs(C-R)+abs(C-L)+abs(R-p[2])+abs(L-p[4])+abs(R-p[1])+abs(L-p[5])
}
kabu3 <- function(par){
kabutogani3(par[1],par[2],par[3],par[4],par[5],par[6])
}
p=ngon(5,axis=T,col='lightblue',lwd=2)
opt=optim(runif(6),kabu3,method='Nelder-Mead')
kabu3(opt$par)
(par=opt$par)
L=par[1]+par[2]*1i
C=par[3]+par[4]*1i
R=par[5]+par[6]*1i
pt(C,'C') ; pt(L,'L') ; pt(R,'R')
seg(p[1],R);seg(p[2],R);seg(p[3],C);seg(p[4],L);seg(p[5],L);seg(L,C);seg(R,C)
0621卵の名無しさん
垢版 |
2018/08/17(金) 09:24:47.60ID:qXmvzl8y
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),paste('p',1:n),cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),col=2,lwd=2,...)
pt <- function(x,y,...) text(Re(x),Im(x), y,...)

kabutogani3 <- function(xl,yl,xc,yc,xr,yr){
L=xl+yl*1i
C=xc+yc*1i
R=xr+yr*1i
abs(p[3]-C)+abs(C-R)+abs(C-L)+abs(R-p[2])+abs(L-p[4])+abs(R-p[1])+abs(L-p[5])
}
kabu3 <- function(par){
kabutogani3(par[1],par[2],par[3],par[4],par[5],par[6])
}
p=ngon(5,axis=T,col='lightblue',lwd=2)
opt=optim(runif(6),kabu3,method='CG')
kabu3(opt$par)
(par=opt$par)
L=par[1]+par[2]*1i
C=par[3]+par[4]*1i
R=par[5]+par[6]*1i
pt(C,'C') ; pt(L,'L') ; pt(R,'R')
seg(p[1],R);seg(p[2],R);seg(p[3],C);seg(p[4],L);seg(p[5],L);seg(L,C);seg(R,C)
0622卵の名無しさん
垢版 |
2018/08/17(金) 17:56:05.86ID:qXmvzl8y
# how many ways of allocating 5 rooms to 6 people without vacancy?

# allocated to 1 room (4 vacant)
a1=choose(5,1)*1^6 ; a1

# allocated to 2 rooms (3 vacant)
a2=choose(5,2)*(2^6-2) ; a2

# allocated to 3 rooms (2 vacant)
a3=choose(5,3)*( 3^6-choose(3,2)*(2^6-2)-3 ) ; a3

# allocated to 4 rooms (1 vacant)
a4=choose(5,4)*( 4^6 - choose(4,3)*(3^6-choose(3,2)*(2^6-2)-3) - choose(4,2)*(2^6-2)-4 ) ; a4

5^6 - a1 - a2 - a3 - a4
0623卵の名無しさん
垢版 |
2018/08/17(金) 18:38:19.94ID:qXmvzl8y
dec2n <- function(num, N = 2, digit = 1){ # decimal to 0,1,..,n-1 vector
if(num <= 0 & digit <= 0) return()
else{
return(append(dec2n(num%/%N, N ,digit-1), num%%N))
}
}

room.allocation <- function(n,r){ # allocate n people to r rooms without vacancy
max=r^n
counter=0
x=0
while(x < max){
if(length(unique(dec2n(x,r,n))) == r) counter = counter+1
x= x + 1
}
return(counter)
}
0624卵の名無しさん
垢版 |
2018/08/18(土) 16:03:38.91ID:Z1UCnKoz
kagamimochi = function(b, h){
r = (b^2 + h^2)/(2*h)
V = (2/3*r^3 - r^2*(r-h) + (r-h)^3/3)*pi
return(c(radius=r,Volume=V))
}
0625卵の名無しさん
垢版 |
2018/08/18(土) 17:38:21.87ID:Z1UCnKoz
kagamimochi = function(b, h){
r = (b^2 + h^2)/(2*h)
if(b > h) V = (2/3*r^3 - r^2*(r-h) + (r-h)^3/3)*pi
else V = (r^2*(h-r) - 1/3*(h-r)^3 + 2*r^3/3)*pi
return(c(radius=r,Volume=V))
}
0626卵の名無しさん
垢版 |
2018/08/18(土) 20:09:40.32ID:TnH4H/8z
draft
# how many ways of allocating 6 rooms to n people without vacancy?

# allocated to 1 room
a1=choose(6,1)*1^n

# allocated to 2 rooms
a2=choose(6,2)*(2^n-2)

# allocated to 3 rooms
a3=choose(6,3)*( 3^n-choose(3,2)*(2^n-2)-3 )

# allocated to 4 rooms
a4=choose(6,4)*( 4^n - choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4 )

# allocated to 5 rooms
a5=choose(6,5)*( 5^n - choose(5,4)*(4^n-choose(4,3)*(3^n-
choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4)-choose (5,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose (5,2)*(2^n -2) -5


6^n - a1 - a2 - a3 - a4 - a5
0627卵の名無しさん
垢版 |
2018/08/18(土) 20:22:13.81ID:Z1UCnKoz
# how many ways of allocating 6 rooms to n people without vacancy?
library(Rmpfr)
six_rooms <- function(x){
if(x[1]<10) n=x else n=mpfr(x,100)
# allocated to 1 room
a1=choose(6,1)*1^n
# allocated to 2 rooms
a2=choose(6,2)*(2^n-2)
# allocated to 3 rooms
a3=choose(6,3)*( 3^n-choose(3,2)*(2^n-2)-3 )
# allocated to 4 rooms
a4=choose(6,4)*( 4^n - choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4 )
# allocated to 5 rooms
a5=choose(6,5)*(5^n-choose(5,4)*(4^n-choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3)
-choose(4,2)*(2^n-2)-4)-choose(5,3)*(3^n-choose(3,2)*(2^n-2)-3)-choose(5,2)*(2^n -2) -5)
6^n - a1 - a2 - a3 - a4 - a5
}
0628卵の名無しさん
垢版 |
2018/08/19(日) 00:27:23.03ID:gaGSkZ47
allocate.rooms <- function(m,n){ # m:rooms n:people
if(m==n) return(factorial(m))
else if(m==1) return(1)
else m*Recall(m,n-1) + m*Recall(m-1,n-1)
}

#include<stdio.h>

long factorial(long n) {
long re = 1;
long k;
for(k=1;k <=n;k++) {re *= k;}
return re;
}

long rooms(int m, int n){
if(m==n) { return factorial(m);}
else if(m==1){ return 1;}
else{
return m * rooms(m,n-1) + m * rooms(m-1,n-1);
}
}

void main( int argc, char *argv[] ){
int m,n;
long ways;
m=atoi(argv[1]);
n=atoi(argv[2]);
ways=rooms(m,n);
printf("%d\n",ways);
}
0629卵の名無しさん
垢版 |
2018/08/19(日) 15:49:26.54ID:gaGSkZ47
楕円体　x^2/a^2 + y^2/b^2 + z^2/c^2 = 1

高さ　h のロケットおっぱいの体積　

RocketPi <= function(a,b,c,h) 2/3*pi*a*b*c - 1/3*pi*b*c*(3*a^2-(a-h)^2)*(a-h)/a
0630卵の名無しさん
垢版 |
2018/08/19(日) 15:56:15.35ID:gaGSkZ47
楕円体　x^2/a^2 + y^2/b^2 + z^2/c^2 = 1

RocketPi <= function(a,β,γ,h)　# a:楕円体の長軸長,　β,γ:ロケットおっぱいの楕円底面の軸長,ｈ：ロケットおっぱいの高さ
b=(a-h)*β/sqrt(a^2-(a-h)^2)
c=(a-h)*γ/sqrt(a^2-(a-h)^2)
2/3*pi*a*b*c - 1/3*pi*b*c*(3*a^2-(a-h)^2)*(a-h)/a
}
0631卵の名無しさん
垢版 |
2018/08/19(日) 17:09:55.54ID:gaGSkZ47
# x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
(1+sqrt(5))/2
f <- function(a=10,b=10,c=10){
xy2z <- function(x,y) c*sqrt(a^2*b^2-a^2*y^2-b^2*x^2)/(a*b)
x=seq(-a,a,le=50)
y=seq(-r*a,r*a,le=50)
z=outer(x,y,xy2z)
contour(x,y,z)
persp(x,y,z,
theta=35, lty=3,col="pink",xlab='x',ylab='y',zlab='z',ylim=c(-2*a,2*a),
ticktype='detailed',shade=0.75,phi=30,ltheta=-10,border=TRUE)
rgl::persp3d(x,y,z,col='pink')
}
f()
f(10,20,25)
0632卵の名無しさん
垢版 |
2018/08/19(日) 23:40:10.33ID:gaGSkZ47
seg <- function(a,b,...){
segments(Re(a),Im(a),Re(b),Im(b),...)}
pt <- function(x,y=NULL,...){
text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)}
# solve α^2/a^2 + (b+β)^2/b^2=1 for a
α=3
β=-1
b=1
# α^2*b^2/(-β*(2*b+β))
(a = α*b/sqrt(-β*(2*b+β)))
(q=b+β)
f= function(x,y) x^2/a^2 + (y-q)^2/b^2 # = 1
x=seq(-5,5,length=100)
y=seq(-5,5,length=100)
z=outer(x,y,f)
contour(x,y,z,level=1,bty='l',xlim=c(-10,10),ylim=c(-2,10),
drawlabels=FALSE,col=sample(colours(),1),axes=FALSE,lwd=2)
pt(α,'●') ; pt(-α,'●') ; pt(β*1i,'●')
seg(-10,10,lty=3) ; seg(-4i,10i,lty=3)
pt(0+q*1i,'.')
eclipse <- function(b) {
# α^2*b^2/(-β*(2*b+β))
(a = α*b/sqrt(-β*(2*b+β)))
(q=b+β)
f= function(x,y) x^2/a^2 + (y-q)^2/b^2 # = 1
x=seq(-10,10,length=100)
y=seq(-10,10,length=100)
z=outer(x,y,f)
contour(x,y,z,level=1,drawlabels=FALSE,add=TRUE,col=sample(colours(),1),lwd=2)
pt(0+q*1i,'.')
}
for(i in 2:10) eclipse(i)
■ このスレッドは過去ログ倉庫に格納されています
ニューススポーツなんでも実況