18/06/30 15:09:55.32 j2CU1Lw0.net
a=750 ; b=250
v=c(rep(1,a),rep(0,b))
f <- function(v){
x=sample(v)
any(cumsum(x)==cumsum(1-x))
}
mean(replicate(1e5,f(v)))
598:卵の名無しさん
18/07/01 10:28:46.39 F3KIhVUK.net
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x)
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}
#
birth(100, 3)
vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
599:卵の名無しさん
18/07/01 11:53:44.57 F3KIhVUK.net
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x))
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}
#
birth(100, 3)
vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
600:卵の名無しさん
18/07/01 12:03:30.64 F3KIhVUK.net
date=1:366
p=c(97/400,rep(1,365))
same.birth <- function(n,lwr=2,upr=1e6){
x=sample(date,n,replace=TRUE,prob=p)
di=max(table(x)
lwr<=di & di<=upr
}
birth <- function(n,lwr=2,upr=1e6,k=1e4){
mean(replicate(k,same.birth(n,lwr,upr)))
}
#
birth(100, 3)
vrb=Vectorize(birth)
x=1:50
y=vrb(x)
plot(x,y,pch=19)
abline(h=0.5,lty=3,col=4)
min(x[whic(y > 0.5)])
601:卵の名無しさん
18/07/01 12:08:34.56 F3KIhVUK.net
インフルエンザの迅速キットは特異度は高いが感度は検査時期によって左右される。
ある診断キットが開発されたとする。
このキットは特異度は９９％と良好であったが、
感度については確かな情報がない。
事前確率分布として一様分布を仮定する。
５０人を無作為抽出してこの診断キットで診断したところ４０人が陽性であった。
この母集団の有病率の期待値と９５％信用区間はいくらか？
またこの診断キットの感度の期待値と９５％信用区間はいくらか
暇つぶしにこれをMCMCを使わずに解く方法を考えていた。
偽陽性率FP=0.01として
陽性確率p=TP*x+(1-x)*FP
尤度が50C40*p^40*(1-p)^10
TPは一様分布なので積分消去して
確率密度関数に比例する関数を作ってarea under the curveで割って確率密度関数化したのち積分して累積密度関数をつくる。この累積密度関数の逆関数を作って95%区間が最短になる区間を計算すれば信頼区間が算出できる。
この結果がstanでのシミュレーションの結果と一致すればよし。
602:卵の名無しさん
18/07/01 22:00:46.73 lSk1MYzX.net
# # choose(n,r) == gamma(n+1) / (gamma(r+1) * gamma(n-r+1))
same.birthday <- function(n) 1-choose(365+97/400,n)*factorial(n)/(365+97/400)^n
plot(x,y,bty='l',xlab='subjects',ylab='probability')
curve(same.birthday(x),add = TRUE)
abline(h=0.5,col=8)
same.birthday(22:23)
603:卵の名無しさん
18/07/02 01:26:48.74 pp47QgIN.net
library(rjags)
N=50
z=40
FP=0.01
shape1=1
shape2=1
dataList=list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
modelstring <- paste0("
model{
theta=TP*x+FP*(1-x)
z ~ dbinom(theta,N)
TP ~ dbeta(shape1,shape2)
x ~ dbeta(shape1,shape2)
}"
)
writeLines( modelstring , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList, quiet=TRUE)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("TP","x","theta"), n.iter=100000 )
js=as.matrix(codaSamples)
head(js)
BEST::plotPost(js[,'TP'],xlab='sensitivity')
BEST::plotPost(js[,'x'],xlab='prevalence')
BEST::plotPost(js[,'theta'],xlab='positive result',showMode = TRUE)
604:卵の名無しさん
18/07/02 01:28:46.92 pp47QgIN.net
N=50
z=40
FP=0.01
shape1=1
shape2=1
data = list(N=N,z=z,FP=FP,shape1=shape1,shape2=shape2)
stanString=paste0('
data{
int N;
int z;
real FP;
real shape1;
real shape2;
}
parameters{
real<lower=0,upper=1> TP;
real<lower=0,upper=1> x; //prevalence
}
transformed parameters{
real<lower=0,upper=1> theta;
theta=TP*x+FP*(1-x);
}
model{
z ~ binomial(N,theta);
TP ~ beta(shape1,shape2); // T[0.5,];
x ~ beta(shape1,shape2);
}
')
605:卵の名無しさん
18/07/02 01:29:10.57 pp47QgIN.net
# model=stan_model(model_code = stanString)
# saveRDS(model,'quick_kit.rds')
model=readRDS('quick_kit.rds')
fit=sampling(model,data=data,iter=10000)
print(fit,digits=3,probs=c(.025,.50,.975))
stan_trace(fit)
# stan_diag(fit)
stan_ac(fit)
stan_dens(fit,separate_chains = TRUE)
stan_hist(fit,fill='skyblue',bins=15,pars=c('x','TP'))
ms=rstan::extract(fit)
BEST::plotPost(ms$TP,showMode = TRUE,xlab='TP')
BEST::plotPost(ms$x,showMode = FALSE,xlab='prevalence',col=sample(colours(),1))
606:卵の名無しさん
18/07/02 21:40:12.88 tMXA02ZR.net
escape.cliff <- function(p=2/3,k=1000){
pos=1
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1)
607:,1,prob=c(p,1-p)) } return(TRUE) } mean(replicate(1e2,escape.cliff()))
608:卵の名無しさん
18/07/02 22:58:14.27 tMXA02ZR.net
escape.cliff <- function(pos=1,p=2/3,k=1000){
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e3,escape.cliff(10,0.5,100)))
#
totter <- function(p=2/3,pos=1,pos0=NULL,k=1e3){
if(is.null(pos0)) pos0=pos
for(i in 1:k) {
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(pos>pos0)
}
mean(replicate(1e3,totter(0.5,10,5 )))
609:卵の名無しさん
18/07/03 07:46:06.00 22JZXDLY.net
exam <- function(p=0.5,hit=0,money=5,k=5){
while(hit<k){
money=money-1
if(money==0) return(FALSE)
shoot=sample(c(1,0),1,prob=c(p,1-p))
if(shoot){
money=money+1
hit=hit +1
}
}
return(TRUE)
}
mean(replicate(1e5, exam()))
610:卵の名無しさん
18/07/03 08:47:34.93 ltuSOSv2.net
ryunen <- function(p=0.6 ,money=10,grade=0,ryu=0){
while(grade<6){
test=rbinom(1,1,p)
if(test) grade=grade+1
else ryu=ryu+1
if(ryu==2){
money=money-5
ryu=ryu-1
}
if(money<=0) return(FALSE)
}
return(grade==6)
}
mean(replicate(1e5,ryunen()))
611:卵の名無しさん
18/07/03 11:46:57.62 22JZXDLY.net
escape.cliff <- function(pos=1,p=2/3,k=1000){
while(pos<k){
if(pos==0) return(FALSE)
pos=pos+sample(c(1,-1),1,prob=c(p,1-p))
}
return(TRUE)
}
mean(replicate(1e5,escape.cliff(1,2/3,3)))
4/7
612:卵の名無しさん
18/07/03 13:41:51.20 22JZXDLY.net
gmbl <- function(money=20,p=18/38){
while(0<money & money <40){
money=money+sample(c(1,-1),1,p=c(p,1-p))
}
return(money==40)
}
mean(replicate(1e3, gmbl()))
613:卵の名無しさん
18/07/03 14:41:00.85 22JZXDLY.net
x=c(rep(1,4),rep(0,48))
mean(replicate(1e5,which(sample(x)==1)[1]))
y=numeric(48)
for(i in 1:48) y[i]=i*choose(48,i-1)*factorial(i-1)*4*factorial(52-i)/factorial(52)
sum(y)
614:卵の名無しさん
18/07/03 14:42:24.05 22JZXDLY.net
(a) A railroad numbers its locomotives in order, 1, 2, . . . , N. One day you see a locomotive and its number is 60.
Guess how many locomotives the company has.
(b) You have looked at 5 locomotives and the largest number observed is 60.
Again guess how many locomotives the company has.
615:卵の名無しさん
18/07/03 16:06:53.54 22JZXDLY.net
n=5
m=60
N=n:100
loco <- function(x){
max(sample(1:x,n))
}
vloco=Vectorize(loco)
loco.sim <- function() N[which(vloco(N)==m)]
locomotives=unlist(replicate(1e4,loco.sim()))
BEST::plotPost(locomotives,breaks=30)
HDInterval::hdi(locomotives)
616:卵の名無しさん
18/07/03 18:26:54.88 22JZXDLY.net
n=5
library(gtools)
y=permutations(n,n,1:n)
f=function(x) sum(x==1:n)
sum(apply(y,1,f))/length(y)
1/n
617:卵の名無しさん
18/07/03 18:39:40.14 22JZXDLY.net
f = function(n){
sum(sample(1:n)==1:n)/n
}
g=function(n)mean(replicate(1e4,f(n)))
vg=Vectorize(g)
x=1:100
plot(x,vg(x))
curve(1/x,add=T)
618:卵の名無しさん
18/07/03 21:26:39.83 22JZXDLY.net
gr=expand.grid(1:9,1:9,1:9)
u=10^(2:0)
f=function(a,b,c,x=1776) sum(c(a,b,c)*u+c(b,c,a)*u+c(c,a,b)*u)-x
# when A>B>C
idx=which(mapply(f,gr[,1],gr[,2],gr[,3])==0 & gr[,1]>gr[,2] & gr[,2]>gr[,3])
gr[idx,]
# otherwise
idx2=which(mapply(f,gr[,1],gr[,2],gr[,3])==0)
length (idx2)
gr[idx2,]
619:卵の名無しさん
18/07/04 00:37:12.00 sDwIwpA3.net
n=5
m=60
N=n:100
loco <- function(x){
max(sample(1:x,n))
}
vloco=Vectorize(loco)
loco.sim <- function() N[which(vloco(N)==m)]
locomotives=unlist(replicate(1e4,loco.sim()))
BEST::plotPost(locomotives,breaks=30)
HDInterval::hdi(locomotives)
n=60:100
pmf=choose(59,4)/choose(n,5) #Pr(max=60|n)
pdf=pmf/sum (pmf)
sum( n*pdf) #E(n)
lines(n,pdf)
620:卵の名無しさん
18/07/05 10:51:42.59 9YHvFei/.net
own <- function(n){
mean(replicate(1e5,sum(sample(1:n)-1:n==0)))
}
own(100)
621:卵の名無しさん
18/07/05 10:54:03.43 9YHvFei/.net
r1 <- function(x){ # rotate by one bead
n=length(x)
y=numeric(n)
y[1]=x[n]
for(i in 1:(n-1)){
y[i+1]=x[i]
}
return(y)
}
rn <- function(x){ # every rotation
n=length(x)
mat=matrix(rep(NA,n^2),ncol=n)
mat[,1]=x
for(i in 1:(n-1)){
mat[,i+1]=r1(mat[,i])
}
return(t(mat))
}
same <- function(x,y){
if(sum(x)!=sum(y)) return(FALSE)
f=function(a,b=y){ # is equal to y
all(a==b)
}
mat=rbind(rn(x),rn(rev(x))) # with symmetric conversion
any(apply(mat,1,f))
}
622:卵の名無しさん
18/07/05 10:55:44.62 9YHvFei/.net
dec2bin <- function(num, digit=0){ # decimal to 0,1 vector
if(num <= 0 && digit <= 0){
return(NULL)
}else{
return(
623:append(Recall(num%/%2,digit-1), num%%2)) } } vd2b=Vectorize(dec2bin) # bracelett <- function(n){ mat=t(vd2b(0:(2^n-1),n)) # make all permutation of beads # head(mat) ret=list() # list of the same bracelett for(i in 1:2^n){ ret[[i]]=which(apply(mat,1,function(z)same(z,mat[i,]))) } # head(ret) ; table(unlist(ret)) del=NULL for(i in 1:2^n) del=append(del,ret[[i]][-1]) 2^n - length(unique(del)) }
624:卵の名無しさん
18/07/05 14:59:01.13 bn3tqJqU.net
> library(gtools)
> swap <- function(n){
+ perm=permutations(n,n,1:n)
+ mean(apply(perm,1,function(x)sum(x==1:n)))
+ }
> swap(9)
[1] 1
> swap.sim <- function(n,k=1e5){
+ mean(replicate(k,sum(sample(1:n)==1:n)))
+ }
> swap.sim(9)
[1] 1.0004
625:卵の名無しさん
18/07/05 20:18:51.37 bn3tqJqU.net
swap2 <- function(n){ # k*nCk*(1/n)^k*(1-1/n)^(n-k) k=0:n
re=numeric(n)
for(i in 0:n){
p=1/n
re[i]=i*choose(n,i)*p^i*(1-p)^(n-i)
}
sum(re)
}
swap2(7)
#
swap3 <- function(n){
re=numeric(n)
for(i in 0:n){
p=1/n
re[i]=i*dbinom(i,n,p)
}
sum(re)
}
swap3(7)
626:卵の名無しさん
18/07/06 09:11:30.20 UEcW6fma.net
f012 <- function(n){
args=list()
length(args)=n-1
args[[1]]=args[[n-1]]=1:2
for(i in 2:(n-2)){
args[[i]]=0:2
}
gr=do.call(expand.grid,args)
gr=as.matrix(gr)
ret=numeric()
for(i in 1:nrow(gr)){
ret[i]=all(diff(gr[i,])!=0)
}
sum(ret)
}
627:卵の名無しさん
18/07/21 16:31:14.15 V1Wm3iKf.net
/* SEND+MORE=MONEYの覆面暗算をC言語で解いてみる。*/
#include<stdio.h> /* SEND+MORE=MONEY */
int compare_int(const void *a, const void *b)
{return *(int*)a - *(int*)b;}
int i,j;
int unique(int num[8]){
qsort(num,8,sizeof(int),compare_int);
for(i=0;i<8;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;
}
main(){
int S,E,N,D,M,O,R,Y;
for(S = 1; S < 10; S++){
for(E = 0; E < 10; E++){
for(N = 0; N < 10; N++){
for(D = 0; D < 10; D++){
for(M = 1; M < 10; M++){
for(O = 0; O < 10; O++){
for(R = 0; R < 10; R++){
for(Y = 0; Y < 10; Y++){
if(S*1000+E*100+N*10+D+M*1000+O*100+R*10+E==M*10000+O*1000+N*100+E*10+Y){
int num[]={S,E,N,D,M,O,R,Y};
if(unique(num)==1){
printf("%d%d%d%d+%d%d%d%d=%d%d%d%d%d\n", S,E,N,D,M,O,R,E,M,O,N,E,Y);
}}}}}}}}}}}
C:\MinGW>gcc sendmore.c -o money
C:\MinGW>money
9567+1085=10652
628:卵の名無しさん
18/07/22 11:46:07.83 tW9s1ZUi.net
#include<stdio.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;}
int unique(int num[]){
int i,j,n=10;
qsort(num,n,sizeof(int),compare_int);
for(i=0;i<n;i++){
for(j=0;j<i;j++){
if(num[j]==num[j+1]){
return 0;
}}}
return 1;}
main(){ /* ド底辺＋私立医＝裏口馬鹿 */
int n=1,A,B,C,D,E,F,G,H,I,J;
for(A = 1; A < 10; A++){
for(B = 0; B < 10; B++){
for(C = 0; C < 10; C++){
for(D = 1; D < 10; D++){
for(E = 0; E < 10; E++){
for(F = 0; F < 10; F++){
for(G = 1; G < 10; G++){
for(H = 0; H < 10; H++){
for(I = 0; I < 10; I++){
for(J = 0; J < 10; J++){
if(A*100+B*10+C +D*100+E*10+F==G*1000+H*100+I*10+J){
int num[]={A,B,C,D,E,F,G,H,I,J};
if(unique(num)==1){
printf("%2d: %d%d%d + %d%d%d = %d%d%d%d\n", n,A,B,C,D,E,F,G,H,I,J);
n++;
}}}}} }}} }}} }}
629:卵の名無しさん
18/07/23 16:11:37.39 B1Q5jsGg.net
# ２つの整数があります。
# それらをたしてできた数は、
# 十の位と一の位の数字が等しい２けたの整数になり、
# それらをかけてできた数は、
# 百の位、十の位、一の位が等しい３けたの整数になりました。
# このような２つの整数の組をあるだけ答えなさい。
gr=expand.grid(1:99,1:99)
f <- function(x,y){ all(x>=y,c((x+y)%%10==(x+y)%/%10,(x*y)%/%100==((x*y)%/%10)%%10,(x*y)%/%100==(x*y)%%10))}
gr[which(mapply(f,gr[,1],gr[,2])==TRUE),]
630:卵の名無しさん
18/07/23 17:34:03.68 B1Q5jsGg.net
# 実行速度遅すぎて実用性なし
# Consider a function which, for a given whole number n,
# returns the number of ones required when writing out all numbers between 0 and n.
# For example, f(13)=6. Notice that f(1)=1.
# What is the next largest number n such that f(n)=n.
ｆ０ ＜－ ｆｕｎｃｔｉｏｎ（ｎ）｛
ｙ＝ａｓ．ｃｈａｒａｃｔｅｒ（ｎ）
ｎｃ＝ｎｃｈａｒ（ｙ）
ｚ＝ＮＵＬＬ
ｆｏｒ（ｉ ｉｎ １：ｎｃ） ｚ［ｉ］ ＜－ ｓｕｂｓｔｒ（ｙ，ｉ，ｉ）＝＝１
ｒｅｔｕｒｎ（ｓｕｍ（ｚ））
｝
ｆ ＜－ ｆｕｎｃｔｉｏｎ（ｎ）｛
ｒｅ＝ｎｕｍｅｒｉｃ（）
ｒｅ［１］＝１
ｆｏｒ（ｉ ｉｎ ２：ｎ）｛
ｒｅ［ｉ］＝ｒｅ［ｉ－１］＋ｆ０（ｉ）
｝
ｒｅｔｕｒｎ（ｒｅ［ｎ］）
｝
ｇ ＜－ ｆｕｎｃｔｉｏｎ（ｎ） ｎ － ｆ（ｎ）
ｉ＝１９９９７９
ｗｈｉｌｅ（ｇ（ｉ）！＝０）｛
ｉ＜－ｉ＋１
｝
ｉ
631:卵の名無しさん
18/07/23 22:11:59.00 B1Q5jsGg.net
6!=(5+1)5!=5*5!+5!
=5*5!+(4+1)*4!
=5*5!+4*4!+4!
=5*5!+4*4!+(3+1)*3!
=5*5!+4*4!+3*3!+3!
=5*5!+4*4!+3*3!+(2+1)*2!
=5*5!+4*4!+3*3!+2*2!+2!
=5*5!+4*4!+3*3!+2*2!+1*1!+1
factorial(n+1)= sum(factoral(n:1)*(n:1))+1
632:卵の名無しさん
18/07/24 09:34:20.21 RZUDxWdB.net
Ｂ，Ｃ，Ｄ，Ｅ，Ｆが０～９の数字（同じ数字であってもよい）で
６！＊Ｂ＋５！＊Ｃ＋４！＊Ｄ＋３！＊Ｅ＋２！＊Ｆ＋１！＊Ｇ＝５５５５
が成立するときＢ＋Ｃ＋Ｄ＋Ｅ＋Ｆ＋Ｇの最小となるＢ～Ｇの組み合わせを求めよ。
n=5555
f <- function(B,C,D,E,F,G) sum(factorial(6:1)*c(B,C,D,E,F,G))-n
max=n%/%factorial(6)
gr=expand.grid(0:max,0:6,0:5,0:4,0:3,0:2)
ret=mapply(f, gr[,1],gr[,2],gr[,3],gr[,4],gr[,5],gr[,6])
bg=apply(gr[which(ret==0),],1,sum)
(min.sum=bg[which.min(bg)])
gr[names(min.sum),]
n%/%factorial(6)
n%%factorial(6)%/%factorial(5)
n%%factorial(6)%%factorial(5)%/%factorial(4)
n%%factorial(6)%%factorial(5)%%factorial(4)%/%factorial(3)
n%%factorial(6)%%factorial(5)%%factorial(4)%%factorial(3)%/%factorial(2)
n%%factorial(6)%%factorial(5)%%factorial(4)%%factorial(3)%%factorial(2)%/%factorial(1)
633:卵の名無しさん
18/07/24 16:06:02.90 Cp/BkN3V.net
#include<stdio.h>
#include<string.h>
int compare_int(const void *a, const void *b){
return *(int*)a - *(int*)b;
}
main( int argc, char *argv[] ){
char N = '0'+ atoi(argv[1]);
int a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t;
for(a = 0; a < 10; a++){
for(b = 0; b < 10; b++){
for(c = 0; c < 10; c++){
for(d = 0; d < 10; d++){
for(e = 0; e < 10; e++){
for(f = 0; f < 10; f++){
for(g = 0; g < 10; g++){
for(h = 0; h < 10; h++){
for(i = 0; i < 10; i++){
for(j = 0; j < 10; j++){
for(k = 0; k < 10; k++){
for(l = 0; l < 10; l++){
634:卵の名無しさん
18/07/24 16:07:14.87 Cp/BkN3V.net
for(m = 0; m < 10; m++){
for(n = 0; n < 10; n++){
for(o = 0; o < 10; o++){
for(p = 0; p < 10; p++){
for(q = 0; q < 10; q++){
for(r = 0; r < 10; r++){
for(s = 0; s < 10; s++){
for(t = 0; t < 10; t++){
int num[]={a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t};
qsort(num,20,sizeof(int),compare_int);
char NUM[20];
int i;
for(i=0;i<20;i++){
NUM[i]='0'+ num[i];
}
char N4[4]={N,N,N,N};
char N5[5]={N,N,N,N,N};
if(strstr(NUM,N4)!=NULL & strstr(NUM,N5)==NULL &
(100*a+10*b+c)*f == 100*g+10*h+i &
(100*a+10*b+c)*e == 100*j+10*k+l &
(100*a+10*b+c)*d == 100*m+10*n+o &
(100*m+10*n+o)*100 + (100*j+10*k+l)*10 + 100*g+10*h+i == p*10000+q*1000+r*100+s*10+t
){
printf("abc def ghi jkl mno pqrst = %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d %d%d%d%d%d\n",
a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t);
}
}}}}} }}}}} }}}}} }}}}}
}
635:卵の名無しさん
18/07/28 06:54:00.42 Y3H/zTDn.net
options(scipen = 10)
fibon <- function(N){
f=numeric(N)
f[1]=1
f[2]=1
if(N>2){for(i in 3:N){
f[i]=f[i-2]+f[i-1]
}}
return(f[N])
}
636:Vectorize(fibon)(1:15) fibon(50) fibon(100)# wrong!
637:卵の名無しさん
18/07/31 22:47:13.76 SHi6FXfD.net
c([],L,L).
c([X|L1],L2,[X|L3]) :- c(L1,L2,L3).
638:卵の名無しさん
18/08/01 09:03:07.74 PtTCNKBk.net
?- L= [1,2,3,4,5,6,7,8,9],c([_,_,_],X,L),c(Y,[_,_,_],X).
L = [1,2,3,4,5,6,7,8,9],
X = [4,5,6,7,8,9],
Y = [4,5,6] ;
No
639:卵の名無しさん
18/08/01 21:24:01.69 nhdpq/Mi.net
last([Last],Last).
last([_|Rest],Last) :- last(Rest,Last).
% last([1,2,3,4,5],Last).
640:卵の名無しさん
18/08/04 20:38:03.88 dITWU8BY.net
日本人の血液型はＡ,Ｏ,Ｂ,ＡＢの比率が４：３：２：１であるという。
それぞれの血液型の人を最低でも各々４、３、２、１人集めるためには必要な人数の期待値はいくらか？
641:卵の名無しさん
18/08/04 21:45:28.37 Tkj7u78K.net
>>597
blood.type <- function(p,need){
n=length(p)
ABO=NULL
enough <- function(x){
pool=numeric(n)
for(i in 1:n) pool[i]=sum(ABO==i)
all(pool >= need)
}
while(!enough(ABO)){
ABO=append(sample(1:n,1,p=p),ABO)
}
return(length(ABO))
}
p=4:1
need=c(4,3,2,1)
re=replicate(1e5,blood.type(p,need))
mean(re)
642:卵の名無しさん
18/08/04 21:52:26.83 Tkj7u78K.net
blood.samples <- function(p=c(4,3,2,1),need=c(1,1,1,1),k=1e5){
mean(replicate(k,blood.type(p,need)))
}
blood.samples()
blood.samples(need=4:1,k=1e4)
643:卵の名無しさん
18/08/05 07:20:33.69 aHWo4FJr.net
ｓｕｂｓｅｔ（［］，［］）．
ｓｕｂｓｅｔ（［Ｆｉｒｓｔ｜Ｒｅｓｔ］，［Ｆｉｒｓｔ｜Ｓｕｂ１］） ：－ ｓｕｂｓｅｔ（Ｒｅｓｔ，Ｓｕｂ１）．
ｓｕｂｓｅｔ（［＿｜Ｒｅｓｔ］，Ｓｕｂ２） ：－ ｓｕｂｓｅｔ（Ｒｅｓｔ，Ｓｕｂ２）．
？－ ｓｕｂｓｅｔ（［ド底辺，特殊，シリツ医大，イカサマ入試，裏口，馬鹿］，＿恥）， ｗｒｉｔｅｌｎ（＿恥），ｆａｉｌ．
644:卵の名無しさん
18/08/05 18:08:20.35 cLIwu37J.net
URLﾘﾝｸ(i.imgur.com)
645:卵の名無しさん
18/08/05 20:30:16.90 aHWo4FJr.net
#URLﾘﾝｸ(i.imgur.com)
N=5
library(gtools)
perm=permutations(n=2,r=N,v=c(1,5), repeats.allowed = TRUE)
move=t(apply(perm,1,cumsum))
p0=0
P=(p0+move)%%6
q0=2
Q=(q0+move)%%6
is.fe <- function(x,y){ # is first encounter after N sec?
all(x[-N]!=y[-N]) & x[N]==y[N]
}
re=NULL
for(i in 1:nrow(P)){
for(j in 1:nrow(Q)){
re=append(re,is.fe(P[i,],Q[j,]))
}
}
(Ans=mean(re)) ; 3^(N-1)/(2^N)^2
cat(Ans*(2^N)^2,'/',(2^N)^2)
646:卵の名無しさん
18/08/06 07:05:55.55 29wsmNi+.net
is.pe <- function(x,y,M) { # encounter M or more than M times within N sec.
sum(x==y) >= M
}
M=2
re=NULL
for(i in 1:nrow(P)){
for(j in 1:nrow(Q)){
re=append(re,is.pe(P[i,],Q[j,],M))
}
}
(Ans=mean(re))
cat(Ans*(2^N)^2,'/',(2^N)^2)
647:卵の名無しさん
18/08/08 19:27:51.20 OkpyQ+1n.net
URLﾘﾝｸ(ailaby.com)
648:卵の名無しさん
18/08/08 19:31:00.52 OkpyQ+1n.net
hanoi <- function(n,from='A',via='B',to='C'){
if(n >= 1){
Recall(n-1,from,to,via)
cat('move',n,'from',from,' to ',to,'\n')
Recall(n-1,via,from,to)
}
}
hanoi(4)
649:卵の名無しさん
18/08/09 05:35:20.13 mD+P4tQf.net
掛け算を再帰関数で定義する。
VAT <- function(n,m=1.08){
if(n==0) 0
else Recall(n-1) + m
}
VAT(250)
650:卵の名無しさん
18/08/09 06:51:20.87 mD+P4tQf.net
総和を再帰関数で書く
sup <- function(v){
if (length(v)==0) 0
else v[1] + Recall(v[-1])
}
sup(1:10)
651:卵の名無しさん
18/08/09 10:57:33.87 JutZ+/A9.net
# Ackermann
> A <- function(m,n){
+ if(m==0) return(n+1)
+ if(n==0) return(Recall(m-1,1))
+ else return(Recall(m-1,Recall(m,n-1)))
+ }
>
> A(2,4)
[1] 11
> A(3,4)
[1] 125
> A(4,1)
Error: evaluation nested too deeply: infinite recursion / options(expressions=)?
>
652:卵の名無しさん
18/08/09 18:49:05.66 JutZ+/A9.net
累積和
# cumsum with for-loop
cumsumL <- function(v){
n=length(v)
re=numeric(n)
re[1]=v[1]
for(i in 1:(n-1)) re[i+1]=re[i]+v[i+1]
re
}
# cumsum with recursive call
cumsumR <- function(v,res=NULL,i=1){
res[1]=v[1]
if(i==length(v)) return(res)
else{
res[i+1] = res[i] + v[i+1]
Recall(v,res,i+1)
}
}
}
653:卵の名無しさん
18/08/09 20:12:53.52 JutZ+/A9.net
# 10進法をN進法でdigit桁表示する
dec2n <- function(num, N = 2, digit = 0){ # decimal to 0,1,..,n-1 vector
if(num==0 & digit==0) return(0)
if(num <= 0 & digit <= 0) return()
else{
return(append(Recall(num%/%N, N ,digit-1), num%%N))
}
}
> dec2n(0)
[1] 0
> dec2n(11,digit=5)
[1] 0 0 1 0 1 1
> dec2n(9,N=5,digit=3)
[1] 0 0 1 4
> dec2n(1000,N=16)
[1] 3 14 8
654:卵の名無しさん
18/08/09 21:28:39.07 JutZ+/A9.net
>>610 (degugged)
dec2n <- function(num, N = 2, digit = 1){ # decimal to 0,1,..,n-1 vector
if(num <= 0 & digit <= 0) return()
else{
return(append(Recall(num%/%N, N ,digit-1), num%%N))
}
}
dec2n(0)
dec2n(11,digit=5)
dec2n(9,N=5,digit=3)
dec2n(1000,N=16)
#
dec2hex <- function(x){ # decimal to hexa
hex=c(0:9,letters[1:6])
n=length(x)
re=numeric(n)
for(i in 1:n){
if(x[i]==0) re[i]='0'
else re[i]=hex[x[i]+1]
}
cat(re,'\n')
}
dec2hex(c(1,0,15))
dec2hex(dec2n(1000,16))
dec2sexa <- function(x) dec2hex(dec2n(x,16))
655:卵の名無しさん
18/08/14 14:27:35.01 Z2jjlChF.net
draft
X <- function(n,red=10,white=90){
rw=red+white
total=martix(0,nrow=n,ncol=2^n)
p=martix(0,nrow=n,ncol=2^n)
total[1,1:2]=c(rw-1,rw+1)
p[1,1:2]=c(white/rw,red/rw)
if(n > 1){
for(i in 1:n){
li=2^i
total[i+1,1:(2*li)]=c(total[i,1:li]-1, total[i,1:li]+1)
p[i+1,1:2*li)]=c(p[i,1:li]*(total[i,1:li]-red)/total[i,1:li], p[i,1:li]*(red)/total[i,1:li])
}
}
return(sum(p[n,]*total[n,]))
}
656:卵の名無しさん
18/08/15 08:08:12.45 SmB+loM1.net
rm(list=ls())
X = function(n,r=10,w=90){ #n:試行数 r: 赤玉数 w:白玉数
rw=r+w # 試行前総玉数
J=rw+n # 総玉数の上限
# s[i,j] i回試行後に総数がj個である確率の行列
s=matrix(0,nrow=n,ncol=J)
s[1,rw-1]=w/rw ; s[1,rw+1]=r/rw # 1回試行後
if(n > 1){
for(i in 2:n){
for(j in r:J){ # jはr未満にはならない
# if(j==1) s[i,j] = s[i-1,j+1]*(j+1-r)/(j+1)
if(j==J) s[i,j] = s[i-1,j-1] * r/(j-1)
else s[i,j] = s[i-1,j-1] * r/(j-1) + s[i-1,j+1]*(j+1-r)/(j+1)
}
}
}
total=sum((r:J)*s[n,r:J])# n回試行後総数の期待値
white=total-r
return(c(total=total,white=white))
}
> vX=Vectorize(X)
> vX(1000:1010)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
total 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5
white 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5
657:卵の名無しさん
18/08/15 12:17:01.96 eDY+Gzxn.net
# シオマネキ
# URLﾘﾝｸ(i.imgur.com)
H=(sqrt(5+2*sqrt(5)))/2
A=0.5+0i
B=A+r
shiomaneki <- function(x,y){
xy=x+y*1i
(abs(xy-A)+abs(xy-B)+x)*2 + 1
}
x=seq(0,Re(B),length=100)
y=seq(0,H,length=100)
z=outer(x,y,shiomaneki)
contour(x,y,z,nlevels=20,lwd=2,bty='l')
x11() ;persp(x,y,z, theta=35, lty=3,col="lightblue",xlab='x',ylab='y',zlab='z',
ticktype='detailed',shade=0.75,phi=30,ltheta=-10,border=TRUE)
rgl::persp3d(x,y,z,col='blue')
sio <- function(xy) shiomaneki(xy[1],xy[2])
optim(c(0.4,0.4),sio,method='Nelder-Mead')
658:卵の名無しさん
18/08/15 12:17:36.49 eDY+Gzxn.net
# カブトガニ
# URLﾘﾝｸ(i.imgur.com)
H=(sqrt(5+2*sqrt(5)))/2
A=0.5+0i
B=A+r
C=0+H*1i
abs(C-B)
kabutogani <- function(y1,y2){
K=0+y1*1i
G=0.5+y2*1i
abs(K-C) + 2*abs(K-G) + 2*abs(G-A) +2*abs(G-B)
}
y1=seq(0,H,length=100)
y2=seq(0,H,length=100)
z=outer(y1,y2,kabutogani)
contour(y1,y2,
659:z,nlevel=100,lwd=2,xlim=c(0.5,1.5),ylim=c(0,1.2),bty='l') kab <- function(y1y2) kabutogani(y1y2[1],y1y2[2]) optim(c(1,0.5),kab, method = 'Nelder-Mead')
660:卵の名無しさん
18/08/15 15:14:03.13 eDY+Gzxn.net
DOP <- function(n,print=FALSE){ # diagonal length of regular plygon of size length 1
q=2*pi/n
r=cos(q)+1i*sin(q)
p=numeric(n+1)
for(i in 1:(n+1)) p[i]=r^(i-1)
D=NULL
if(n>3){
for(j in 3:(2+ceiling(n-3)/2)){
D= append(D,abs(p[1]-p[j]))
}
}
if(print){
plot(p,type='l',bty='l',axes=FALSE,ann=FALSE,lwd=2)
for(i in 3:(n-1)) segments(1,0,Re(p[i]),Im(p[i]),col='gray')
}
return(D)
}
DOP(17,p=T)
661:卵の名無しさん
18/08/15 22:04:58.96 SmB+loM1.net
DOP <- function(n,print=FALSE){ # diagonal length of regular plygon of size length 1
q=2*pi/n
r=cos(q)+1i*sin(q)
p=numeric(n+1)
for(i in 1:(n+1)) p[i]=r^(i-1)
D=NULL
if(n>3){
for(j in 3:(2+ceiling(n-3)/2)){
D= append(D,abs(p[1]-p[j])/abs(p[1]-p[2]))
}
}
if(print){
plot(p,type='l',bty='l',axes=FALSE,ann=FALSE,lwd=2)
for(i in 3:(n-1)) segments(1,0,Re(p[i]),Im(p[i]),col='gray')
}
return(D)
}
DOP(17,p=T)
662:卵の名無しさん
18/08/15 23:51:51.61 eDY+Gzxn.net
> (sin(2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - (cos(4*2*pi/7)*sin(2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) - sin(4*2*pi/7)/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2) + (cos(2*pi/7)* sin(4*2*pi/7))/((1 - cos(2*pi/7))^2 + sin(2*pi/7)^2))
[1] 2.190643
663:卵の名無しさん
18/08/16 22:05:33.14 8kKe3yXf.net
rm(list=ls())
graphics.off()
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,lwd=1,asp=1)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),1:n,cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
664:卵の名無しさん
18/08/17 09:11:30.66 qXmvzl8y.net
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),paste('p',1:n),cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),col=2,lwd=2,...)
pt <- function(x,y,...) text(Re(x),Im(x), y,...)
kabutogani3 <- function(xl,yl,xc,yc,xr,yr){
L=xl+yl*1i
C=xc+yc*1i
R=xr+yr*1i
abs(p[3]-C)+abs(C-R)+abs(C-L)+abs(R-p[2])+abs(L-p[4])+abs(R-p[1])+abs(L-p[5])
}
kabu3 <- function(par){
kabutogani3(par[1],par[2],par[3],par[4],par[5],par[6])
}
p=ngon(5,axis=T,col='lightblue',lwd=2)
opt=optim(runif(6),kabu3,method='Nelder-Mead')
kabu3(opt$par)
(par=opt$par)
L=par[1]+par[2]*1i
C=par[3]+par[4]*1i
R=par[5]+par[6]*1i
pt(C,'C') ; pt(L,'L') ; pt(R,'R')
seg(p[1],R);seg(p[2],R);seg(p[3],C);seg(p[4],L);seg(p[5],L);seg(L,C);seg(R,C)
665:卵の名無しさん
18/08/17 09:24:47.60 qXmvzl8y.net
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){
r=exp(2*pi/n*1i)
p=complex(n)
for(i in 1:(n+1)) p[i]= (1-r^i)/(1-r)
plot(p,bty='l',type='l',axes=axis, ann=FALSE,asp=1,...)
points(1/(1-r),pch='.')
if(digit) text(Re(p),Im(p),paste('p',1:n),cex=cex)
if(axis){axis(1) ; axis(2)}
invisible(p)
}
seg <- function(a,b,...) segments(Re(a),Im(a),Re(b),Im(b),col=2,lwd=2,...)
pt <- function(x,y,...) text(Re(x),Im(x), y,...)
kabutogani3 <- function(xl,yl,xc,yc,xr,yr){
L=xl+yl*1i
C=xc+yc*1i
R=xr+yr*1i
abs(p[3]-C)+abs(C-R)+abs(C-L)+abs(R-p[2])+abs(L-p[4])+abs(R-p[1])+abs(L-p[5])
}
kabu3 <- function(par){
kabutogani3(par[1],par[2],par[3],par[4],par[5],par[6])
}
p=ngon(5,axis=T,col='lightblue',lwd=2)
opt=optim(runif(6),kabu3,method='CG')
kabu3(opt$par)
(par=opt$par)
L=par[1]+par[2]*1i
C=par[3]+par[4]*1i
R=par[5]+par[6]*1i
pt(C,'C') ; pt(L,'L') ; pt(R,'R')
seg(p[1],R);seg(p[2],R);seg(p[3],C);seg(p[4],L);seg(p[5],L);seg(L,C);seg(R,C)
666:卵の名無しさん
18/08/17 17:56:05.86 qXmvzl8y.net
# how many ways of allocating 5 rooms to 6 people without vacancy?
# allocated to 1 room (4 vacant)
a1=choose(5,1)*1^6 ; a1
# allocated to 2 rooms (3 vacant)
a2=choose(5,2)*(2^6-2) ; a2
# allocated to 3 rooms (2 vacant)
a3=choose(5,3)*( 3^6-choose(3,2)*(2^6-2)-3 ) ; a3
# allocated to 4 rooms (1 vacant)
a4=choose(5,4)*( 4^6 - choose(4,3)*(3^6-choose(3,2)*(2^6-2)-3) - choose(4,2)*(2^6-2)-4 ) ; a4
5^6 - a1 - a2 - a3 - a4
667:卵の名無しさん
18/08/17 18:38:19.94 qXmvzl8y.net
dec2n <- function(num, N = 2, digit = 1){ # decimal to 0,1,..,n-1 vector
if(num <= 0 & digit <= 0) return()
else{
return(append(dec2n(num%/%N, N ,digit-1), num%%N))
}
}
room.allocation <- function(n,r){ # allocate n people to r rooms without vacancy
max=r^n
counter=0
x=0
while(x < max){
if(length(unique(dec2n(x,r,n))) == r) counter = counter+1
x= x + 1
}
return(counter)
}
668:卵の名無しさん
18/08/18 16:03:38.91 Z1UCnKoz.net
kagamimochi = function(b, h){
r = (b^2 + h^2)/(2*h)
V = (2/3*r^3 - r^2*(r-h) + (r-h)^3/3)*pi
return(c(radius=r,Volume=V))
}
669:卵の名無しさん
18/08/18 17:38:21.87 Z1UCnKoz.net
kagamimochi = function(b, h){
r = (b^2 + h^2)/(2*h)
if(b > h) V = (2/3*r^3 - r^2*(r-h) + (r-h)^3/3)*pi
else V = (r^2*(h-r) - 1/3*(h-r)^3 + 2*r^3/3)*pi
return(c(radius=r,Volume=V))
}
670:卵の名無しさん
18/08/18 20:09:40.32 TnH4H/8z.net
draft
# how many ways of allocating 6 rooms to n people without vacancy?
# allocated to 1 room
a1=choose(6,1)*1^n
# allocated to 2 rooms
a2=choose(6,2)*(2^n-2)
# allocated to 3 rooms
a3=choose(6,3)*( 3^n-choose(3,2)*(2^n-2)-3 )
# allocated to 4 rooms
a4=choose(6,4)*( 4^n - choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4 )
# allocated to 5 rooms
a5=choose(6,5)*( 5^n - choose(5,4)*(4^n-choose(4,3)*(3^n-
choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4)-choose (5,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose (5,2)*(2^n -2) -5
6^n - a1 - a2 - a3 - a4 - a5
671:卵の名無しさん
18/08/18 20:22:13.81 Z1UCnKoz.net
# how many ways of allocating 6 rooms to n people without vacancy?
library(Rmpfr)
six_rooms <- function(x){
if(x[1]<10) n=x else n=mpfr(x,100)
# allocated to 1 room
a1=choose(6,1)*1^n
# allocated to 2 rooms
a2=choose(6,2)*(2^n-2)
# allocated to 3 rooms
a3=choose(6,3)*( 3^n-choose(3,2)*(2^n-2)-3 )
# allocated to 4 rooms
a4=choose(6,4)*( 4^n - choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3) - choose(4,2)*(2^n-2)-4 )
# allocated to 5 rooms
a5=choose(6,5)*(5^n-choose(5,4)*(4^n-choose(4,3)*(3^n-choose(3,2)*(2^n-2)-3)
-choose(4,2)*(2^n-2)-4)-choose(5,3)*(3^n-choose(3,2)*(2^n-2)-3)-choose(5,2)*(2^n -2) -5)
6^n - a1 - a2 - a3 - a4 - a5
}
672:卵の名無しさん
18/08/19 00:27:23.03 gaGSkZ47.net
allocate.rooms <- function(m,n){ # m:rooms n:people
if(m==n) return(factorial(m))
else if(m==1) return(1)
else m*Recall(m,n-1) + m*Recall(m-1,n-1)
}
#include<stdio.h>
long factorial(long n) {
long re = 1;
long k;
for(k=1;k <=n;k++) {re *= k;}
return re;
}
long rooms(int m, int n){
if(m==n) { return factorial(m);}
else if(m==1){ return 1;}
else{
return m * rooms(m,n-1) + m * rooms(m-1,n-1);
}
}
void main( int argc, char *argv[] ){
int m,n;
long ways;
m=atoi(argv[1]);
n=atoi(argv[2]);
ways=rooms(m,n);
printf("%d\n",ways);
}
673:卵の名無しさん
18/08/19 15:49:26.54 gaGSkZ47.net
楕円体 x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
高さ h のロケットおっぱいの体積
RocketPi <= function(a,b,c,h) 2/3*pi*a*b*c - 1/3*pi*b*c*(3*a^2-(a-h)^2)*(a-h)/a
674:卵の名無しさん
18/08/19 15:56:15.35 gaGSkZ47.net
楕円体 x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
RocketPi <= function(a,β,γ,h) # a:楕円体の長軸長, β,γ:ロケットおっぱいの楕円底面の軸長,ｈ：ロケットおっぱいの高さ
b=(a-h)*β/sqrt(a^2-(a-h)^2)
c=(a-h)*γ/sqrt(a^2-(a-h)^2)
2/3*pi*a*b*c - 1/3*pi*b*c*(3*a^2-(a-h)^2)*(a-h)/a
}
675:卵の名無しさん
18/08/19 17:09:55.54 gaGSkZ47.net
# x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
(1+sqrt(5))/2
f <- function(a=10,b=10,c=10){
xy2z <- function(x,y) c*sqrt(a^2*b^2-a^2*y^2-b^2*x^2)/(a*b)
x=seq(-a,a,le=50)
y=seq(-r*a,r*a,le=50)
z=outer(x,y,xy2z)
contour(x,y,z)
persp(x,y,z,
theta=35, lty=3,col="pink",xlab='x',ylab='y',zlab='z',ylim=c(-2*a,2*a),
ticktype='detailed',shade=0.75,phi=30,ltheta=-10,border=TRUE)
rgl::persp3d(x,y,z,col='pink')
}
f()
f(10,20,25)
676:卵の名無しさん
18/08/19 23:40:10.33 gaGSkZ47.net
seg <- function(a,b,...){
segments(Re(a),Im(a),Re(b),Im(b),...)}
pt <- function(x,y=NULL,...){
text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)}
# solve α^2/a^2 + (b+β)^2/b^2=1 for a
α=3
β=-1
b=1
# α^2*b^2/(-β*(2*b+β))
(a = α*b/sqrt(-β*(2*b+β)))
(q=b+β)
f= function(x,y) x^2/a^2 + (y-q)^2/b^2 # = 1
x=seq(-5,5,length=100)
y=seq(-5,5,length=100)
z=outer(x,y,f)
contour(x,y,z,level=1,bty='l',xlim=c(-10,10),ylim=c(-2,10),
drawlabels=FALSE,col=sample(colours(),1),axes=FALSE,lwd=2)
pt(α,'●') ; pt(-α,'●') ; pt(β*1i,'●')
seg(-10,10,lty=3) ; seg(-4i,10i,lty=3)
pt(0+q*1i,'.')
eclipse <- function(b) {
# α^2*b^2/(-β*(2*b+β))
(a = α*b/sqrt(-β*(2*b+β)))
(q=b+β)
f= function(x,y) x^2/a^2 + (y-q)^2/b^2 # = 1
x=seq(-10,10,length=100)
y=seq(-10,10,length=100)
z=outer(x,y,f)
contour(x,y,z,level=1,drawlabels=FALSE,add=TRUE,col=sample(colours(),1),lwd=2)
pt(0+q*1i,'.')
}
for(i in 2:10) eclipse(i)
677:卵の名無しさん
18/08/20 13:20:40.68 qQ1lYpXa.net
solve α=a-h, α^2/a^2+β^2/b^2=1,α^2/a^2+γ^2/c^2=1,(a-h+H)^2/a^2 + (β/2)^2/b^2 = 1, for a
678:卵の名無しさん
18/08/21 02:23:50.71 AVDxqhtS.net
rm(list=ls())
D=1
H=0
n=3
cards=c(rep(D,13),rep(H,39))
n.DH=length(cards)
n.D=sum(cards)
sim <- function(){
index_of_inbox=sample(1:n.DH,1)
inbox=cards[index_of_inbox]
outbox=cards[-index_of_inbox] # cards out of box
drawn=sample(outbox,n) # 2 cards drawn from outbox
c(inbox=inbox,drawn=drawn)
}
rate_sim <- function(k){
re=replicate(k,sim())
sum(apply(re,2,function(x) sum(x))==(n+1))/sum(apply(re,2,function(x) sum(x[-1]))==n)
}
re=replicate(100,rate_sim(1000))
summary(re)
#
f <- function(n,D=13,H=52-13){
T=D+H
(D/T * choose(D-1,n)/choose(T-1,n)) /((D/T * choose(D-1,n)/choose(T-1,n)) + H/T * choose(D,n)/choose(T-1,n))
}
n=0:12
f(n)
plot(n,f(n),pch=19,bty='l')
679:卵の名無しさん
18/08/21 12:11:49.73 AVDxqhtS.net
rm(list=ls()) ; graphics.off()
n=5
ngon <- function(n,digit=TRUE,axis=FALSE,cex=1,...){ # draw n-polygon
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) # return vertex complex}
seg <- function(a,b,...){# draw segment of complex a to complex
680:b segments(Re(a),Im(a),Re(b),Im(b),col=2,...)} pt <- function(x,y=NULL,...){ # draw text y at complex x text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)} poly_demo2 <- function(x=runif(n),y=runif(n)){ #各頂点から最短の分岐点を選ぶ x:Re(z), y:Im(z) p=ngon(n,axis=T,col='skyblue') Q=complex(n) re1=re2=0 for(i in 1:n){ # plot Q Q[i]=x[i]+y[i]*1i pt(Q[i],paste('q',i))} for(i in 1:n){ # draw shortest seg for each vertex j=which.min(abs(p[i]-Q)) # Q[j]=nearest node (j: smallest index) seg(p[i],Q[j]) re1=re1+abs(p[i]-Q[j]) # total length of p-Q} for(i in 1:(n-1)){ seg(Q[i],Q[i+1]) re2=re2+abs(Q[i]-Q[i+1]) # total length of Q} return(sum(re1)+sum(re2))}
681:卵の名無しさん
18/08/21 12:13:21.67 AVDxqhtS.net
p=ngon(n,axis=T,col='skyblue')
poly2 <- function(par){ # par=c(Re(z),Im(z))
x=par[1:n]
y=par[(n+1):(2*n)]
Q=complex(n)
re1=re2=0
for(i in 1:n){ # par to complex number
Q[i]=x[i]+y[i]*1i
# pt(Q[i],paste('q',i)) }
for(i in 1:n){
j=which.min(abs(p[i]-Q)) # Q[j]=nearest node (j: smallest index)
# seg(p[i],Q[j])
re1=re1+abs(p[i]-Q[j]) # total length of p-Q}
for(i in 1:(n-1)){
# seg(Q[i],Q[i+1])
re2=re2+abs(Q[i]-Q[i+1]) # total length of Q}
return(sum(re1)+sum(re2))}
par(mfrow=c(2,2))
f <- function(){
(opt=optim(runif(2*n),poly2,method = 'BFGS'))
par=opt$par
poly_demo2(par[1:n],par[(n+1):(2*n)])
}
replicate(4,f())
682:卵の名無しさん
18/08/21 13:27:57.74 AVDxqhtS.net
# n 角形で m 点分岐を最適化(頂点から最短の分岐点を選ぶ)
n=5
m=1
poly_demo3 <- function(x=runif(m),y=runif(m)){ # x:Re(z), y:Im(z)
p=ngon(n,axis=T,col='skyblue')
Q=complex(m)
re1=re2=0
for(i in 1:m){ # plot Q
Q[i]=x[i]+y[i]*1i
pt(Q[i],paste('q',i))
}
for(i in 1:n){ # draw shortest seg for each vertex
j=which.min(abs(p[i]-Q)) # Q[j]=nearest node (j: smallest index)
seg(p[i],Q[j])
re1=re1+abs(p[i]-Q[j]) # total length of p-Q
}
if(m>1){
for(i in 1:(m-1)){
seg(Q[i],Q[i+1])
re2=re2+abs(Q[i]-Q[i+1]) # total length of Q
}}
return(sum(re1)+sum(re2))
}
poly_demo3()
683:卵の名無しさん
18/08/21 13:28:24.29 AVDxqhtS.net
p=ngon(n,print=FALSE)
poly3 <- function(par){ # par=c(Re(z),Im(z))
x=par[1:m]
y=par[(m+1):(2*m)]
Q=complex(m)
re1=re2=0
for(i in 1:m){ # par to complex number
Q[i]=x[i]+y[i]*1i
}
for(i in 1:n){
j=which.min(abs(p[i]-Q)) # Q[j]=nearest node (j: smallest index)
re1=re1+abs(p[i]-Q[j]) # total length of p-Q
}
if(m >1){
for(i in 1:(m-1)){
re2=re2+abs(Q[i]-Q[i+1]) # total length of Q
}}
return(sum(re1)+sum(re2))
}
poly3(runif(6))
par(mfrow=c(2,2))
f <- function(){
(opt=optim(rnorm(2*m),poly3,method = 'CG'))
par=opt$par
poly_demo3(par[1:m],par[(m+1):(2*m)])
}
replicate(4,f())
684:卵の名無しさん
18/08/21 18:03:02.00 AVDxqhtS.net
> str(re)
List of 2
$ :List of 2
..$ value: num 4.02
..$ par : num [1:6] 0.9432 0.6545 0.2652 0.0782 0.4755 ...
$ :List of 2
..$ value: num 3.89
..$ par : num [1:6] 0.8549 0.2805 0.0456 0.3258 0.3861 ...
> re[[which.min( unlist( lapply(re,'[',1) ) )]]
$`value`
[1] 3.891156871891314
$par
[1] 0.85487391200975060 0.28046877582694385 0.04560929326716642
[4] 0.32577844877585493 0.38608966889657914 0.91371285778407096
685:卵の名無しさん
18/08/22 11:10:22.52 rCdgdwPe.net
> dia=1
> heart=0
> n=98
> cards=c(rep(dia,99),rep(heart,1))
> n.DH=length(cards)
> n.D=length(dia)
> sim <- function(){
+ index_of_inbox=sample(1:n.DH,1)
+ inbox=cards[index_of_inbox]
+ outbox=cards[-index_of_inbox] # cards out of box
+ drawn=sample(outbox,n) # n cards drawn from outbox
+ c(inbox=inbox,drawn=drawn)
+ }
> rate_sim <- function(k){
+ re=replicate(k,sim()) # inbox=D&drawn=D / drawn=D
+ all_dia=sum(apply(re,2,function(x) sum(x))==(n+1))
+ drawn_dia=sum(apply(re,2,function(x) sum(x[-1]))==n)
+ c(all_dia/drawn_dia, drawn_dia/k)
+ }
> re=replicate(100,rate_sim(10000))
> summary(re[1,])
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.4235 0.4816 0.5012 0.5001 0.5230 0.5833
> summary(re[2,])
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0170 0.0190 0.0199 0.0199 0.0207 0.0237
686:卵の名無しさん
18/08/22 13:49:25.38 NbYiL2GI.net
GCD <- function(a,b){ # Euclidean mutual division method
r = a%%b # a=bq+r -> a%%b=b%%r
while(r){
a = b
b = r
r = a%%b
}
return(b)
}
reduce_fraction <- function(x,y){
a=x
b=y
r = a%%b # a=bq+r -> a%%b=b%%r
while(r){
a = b
b = r
r = a%%b
}
gcd=b
cat(x/gcd,'/',y/gcd,'\n')
invisible(c(x/gcd,y/gcd))
}
# 2860 / 1082900
reduce_fraction(2860,1082900)
687:卵の名無しさん
18/08/23 16:08:47.89 MpxHuSF8.net
options(digits=16)
a=355
b=113
c=a/b # 3.141592920353983 > pi
d=a/b-pi # 2.667641894049666e-07 > 0
f <- function(i){
j=0 ; x=c ; y=d
while(x < c){
if(y < d) return(j)
x=(a+j)/(b+i)
y=abs(x-pi)
j=j+1
}
return(NULL)
}
i=1
while(is.null(f(i))){
if( !is.null(f(i)) ) return(
688:c(i,f(i))) i=i+1 }
689:卵の名無しさん
18/08/23 17:43:39.29 MpxHuSF8.net
options(digits=9)
a=355
b=113
(U=a/b) # 3.141592 92
(L=2*pi-a/b) # 3.141592 39
f <- function(n){
m=1
x=m/n
while(x<U){
if(L<x) return(m-1)
x=m/n
m=m+1
}
return(NULL)
}
n=1
while(is.null(f(n))){
if( !is.null(f(n)) ) break
n=n+1
}
cat(f(n),'/',n)
690:卵の名無しさん
18/08/23 18:36:03.33 MpxHuSF8.net
LU2fra <- function(L,U){
f <- function(n){
m=0
x=m/n
while(x<U){
if(L<x) return(m-1)
x=m/n
m=m+1
}
return(NULL)
}
n=1
while(is.null(f(n))){
if( !is.null(f(n)) ) break
n=n+1
}
cat(f(n),'/',n,'\n')
f(n)/n
}
dec2fra <- function(digit,precision=0.01){
L=digit*(1-precision)
U=digit*(1+precision)
LU2fra(L,U)
}
691:卵の名無しさん
18/08/25 17:23:40.06 bgFVf9KM.net
Twelve coins are numbered 123456789abc in hex.Balancing coins as follows, 0 denotes logically concluded regular coin. We can use logically concluded regular coins.
(1) 1234=5678 : 0000 0000 9abc
(2) 9ab=000 : 0000 0000 000c c cannot be regular.
(3) c<0 : 0000 0000 000c light
(3) c>0 : 0000 0000 000c heavy
(2) 9ab<000 : 0000 0000 9ab0
(3) 9=a : 0000 0000 00b0 light
(3) 9<a : 0000 0000 9000 light
(3) 9>a : 0000 0000 0a00 light
(2) 9ab>000 : 0000 0000 9ab0
(3) 9=a : 0000 0000 00b0 heavy
(3) 9<a : 0000 0000 0a00 heavy
(3) 9>a : 0000 0000 9000 heavy
(1) 1234<5678 : 1234 5678 0000
Take two light coins and one heavy coin on each dish of the balance
(2) 125=346 : 0000 0078 0000
(3) 7=0 : 0000 0008 0000 heavy
(3) 7>0 : 0000 0070 0000 heavy 7 cannot be in light group
(2) 125<346 : 1200 0600 0000 3,4,and 5 cannot be in both light and heavy group
Balance two coins in light group in (1) and (2),i.e. coin 1 and coin 2
(3) 1=2 : 0000 0600 0000 heavy
(3) 1<2 : 1000 0000 0000 light
(3) 1>2 : 0200 0000 0000 light
(2) 125>346 : 0034 5000 0000 1,2,and 6 cannot ben in both light and heavy group
(3) 3=4 : 0000 5000 0000 heavy
(3) 3<4 : 0030 0000 0000 light
(3) 3>4 : 0004 0000 0000 lihgt
692:卵の名無しさん
18/08/25 17:24:44.77 bgFVf9KM.net
(1) 1234>5678 : 1234 5678 0000
Take two heavy coins and one light coin on each dish of the balance
(2) 125=346 : 0000 0078 0000
(3) 7=0 : 0000 0008 0000 light
(3) 7>0 : 0000 0070 0000 light 7 cannot be in heavy group
(2) 125>346 : 1200 0600 0000 3,4,and 5 cannot be in both heavy and light group
Balance two coins in heavy group in (1) and (2),i.e. coin 1 and coin 2
(3) 1=2 : 0000 0600 0000 light
(3) 1<2 : 1000 0000 0000 heavy
(3) 1>2 : 0200 0000 0000 heavy
(2) 125<346 : 0034 5000 0000 1,2,and 6 cannot ben in both heavy and light group
Balance two coins in heavy group in (1) and (2),i.e. coin 3 and coin 4
(3) 3=4 : 0000 5000 0000 light
(3) 3<4 : 0030 0000 0000 heavy
(3) 3>4 : 0004 0000 0000 lihgt
693:卵の名無しさん
18/08/25 17:56:54.86 LIvDy2Gd.net
6789>abcd
6789<abcd
same as the case with12 coins
(1)6789=abcd 12345 iregular
(2)12=30 45 irregular
(3)4=0 5 irregular
(3)4<0 4 light
(3)4>0 4 heavy
(2)12<30
(3)1=2 3 heavy
(3)1<2 1 light
(3)1>2 2 light
(2)12>30
(3)1=2 3 light
(3)1<2 2 heavy
(3)1>2 1heavy
694:卵の名無しさん
18/08/25 22:00:35.12 LIvDy2Gd.net
deci2frac <- function(x,precision=1e-5){
ae=abs(x*precision)
#accepted error
d=Inf
k=1
while(d>ae){
y=(ifelse(abs(floor(k*x)/k-x) < abs(floor(k*x+1)/k-x),floor(k*x),floor(k*x+1)))
z=y/k
d=abs(x-z)
k=k+1
}
cat(y,'/',k-1,' ',z,'\n')
invisible(c(y,k-1,z))
}
> deci2frac(0.333,precision = 0.01)
1 / 3 0.3333333333333333>
> deci2frac(pi,precision = 10^-7)
355 / 113 3.141592920353983>
> deci2frac(sqrt(2),p=1e-10)
114243 / 80782 1.414213562427273>
> deci2frac(exp(1),1e-8)
20504 / 7543 2.718281850722524>
> exp(1)
[1] 2.718281828459045
>
695:卵の名無しさん
18/08/26 22:49:12.10 o/OlEQv4.net
#include<stdio.h>
void hanoi(int n,char *a,char *b,char *c){
if(n>=1){
hanoi(n-1,a,c,b);
printf("%s から %s へ移す\n",a,c);
/* printf("%d を %s から %s へ移す\n",n,a,c); */
hanoi(n-1,b,a,c);
}
}
void main(int argc,char *argv[]){
int i;
i=atoi(argv[1]);
hanoi(i,argv[2],argv[3],argv[4]);
}
696:卵の名無しさん
18/08/26 23:31:28.09 o/OlEQv4.net
hanoi <- function(n,a='A',b='B',c='C'){
if(n > 0){
Recall(n-1,a,c,b)
cat('move',n,'from',a,' to ',c,'\n')
Recall(n-1,b,a,c)
}
}
options(digits=16)
(2^64-1)/365.2425/24/60/60
library(Rmpfr)
mpfr((2^64-1)/365.2425/24/60/60,1024)
mpfr((2^100-1)/365.2425/24/60/60,1024)
697:卵の名無しさん
18/08/27 22:12:09.93 dp8t1y8e.net
gcd <- function(a,b){
r=a%%b
if(!r) return(b)
gcd(b,r)
}
gcd(18,48)
698:卵の名無しさん
18/08/28 13:27:10.91 tep78gti.net
let rec fact n = if n = 0 then 1 else n * fact (n - 1);;
let rec rms(m,n) = if m=n then fact(m) else if m=1 then 1 else m*rms(m,n-1)+ m*rms(m-1,n-1);;
699:卵の名無しさん
18/08/28 13:41:44.63 tep78gti.net
let rec fact n = if n = 0 then 1 else n * fact (n - 1);;
let rec rms(m,n) = if m=n then fact(m) else if m=1 then 1 else m*rms(m,n-1)+ m*rms(m-1,n-1);;
let ways = rms(int_of_string Sys.argv.(1),int_of_string Sys.argv.(2));;
print_int ways;;
700:卵の名無しさん
18/08/28 16:09:37.85 tep78gti.net
(* nCr = nCr-1 * (n - r + 1) / r *)
let rec nCr(n,r) = if r=0 then 1 else nCr(n,r-1)*(n-r+1)/r;;
(* nCr = n-1Cr-1 + n-1Cr *)
let rec nCr(n,r) = if(n=r || r=0) then 1 else if (n=0 || r> n) then 0 else nCr(n-1,r-1)+ nCr(n-1,r);;
701:卵の名無しさん
18/08/28 16:58:59.85 tep78gti.net
let rec hanoi(n,a,b,c) = if n=1 then print_string("move "^string_of_int(n)^" from "^a^" to "^c^"\n")
else begin
hanoi(n-1,a,c,b);
print_string("move "^string_of_int(n)^" from "^a^" to "^c^"\n");
hanoi(n-1,b,a,c);
end;;
let() = hanoi(int_of_string Sys.argv.(1),Sys.argv.(2),Sys.argv.(3),Sys.argv.(4));;
(* Note ; not ;; between begin and end, begin and end can be ( and ).
to compile
ocamlc.opt -o hanoi.exe hanoi.ml
to execute
hanoi 5 A B C
*)
702:卵の名無しさん
18/08/28 20:20:56.13 ajSYk/Fw.net
let rec facti (n, a) = if n = 0 then a else facti (n - 1, a * n);;
let fact m = facti(m,1);;
fact 10;;
let fact n = let rec facti (n, a) = if n = 0 then a else facti (n - 1, a * n) in facti(n,1);;
fact 10;;
703:卵の名無しさん
18/08/28 20:32:41.05 ajSYk/Fw.net
let rec fibo (n, a1, a2) = if n = 0 then a1 else fibo (n - 1, a1 + a2, a1);;
let fibonacci m = fibo(m,1,0);;
fibonacci 10;;
let fibonacci n = let rec fibo (n, a1, a2) = if n = 0 then a1 else fibo (n - 1, a1 + a2, a1) in fibo (n, 1, 0) ;;
fibonacci 10;;
704:卵の名無しさん
18/08/30 19:17:52.91 9jLRmK1l.net
Ent <- function(P,Q){ # Expected number of trials P:winning Q:losing}
N=P+Q
re=numeric(Q+1)
for(i in 1:(Q+1)) re[i]=i*choose(Q,i-1)/choose(N,i-1)*P/(N-i+1)
return(sum(re))
# Σ[1,Q+1] i * Q!/Q-i+1! * N-i+1!/N! * P/N-i+1!
}
ent <- function(N,p){
P=p*N
Q=(1-p)*N
Ent(P,Q)
}
705:卵の名無しさん
18/08/30 19:18:16.71 9jLRmK1l.net
lottery_sim <- function(N,p,q,hit=1){ # N:lots, p:probability of winning
x=c(rep(1,N*p),rep(0,N*(1-p))) # q:price perlot, hit: hits required to win
win=0
try=0
while(win<hit){
n=length(x)
j=sample(n,1)
win = win + as.numeric(x[j]==1)
x=x[-j]
try=try+1
}
return(try*q)
}
706:卵の名無しさん
18/08/30 20:36:05.34 WRrKmlGk.net
draft
i=5
0011 Q/N*Q-1/N-1*P/N-2*P-1/N-3 ＊ P-2/N-4
0101 Q/N*P/N-1*Q-1/N-2*P-1/N-3
0110
1001
1010
1100
N(N-1)(N-2)(N-3)(N-4)=choose(N,i)*factorial(i)
nPr=function(n,r) choose(n,r)*factorial(r)
i*choose(i-1,hit-1)*nPr(Q,i-hit)*nPr(P,hit)/nPr(N,i)
707:卵の名無しさん
18/08/30 21:07:37.62 WRrKmlGk.net
## Expected Number of Trials 2
# P:atari Q:hazure hit:atari needed to win
Ent2 <- function(P,Q,hit=1){
N=P+Q
nPr=function(n,r) choose(n,r)*factorial(r)
re=numeric(hit+Q)
for(i in hit:(hit+Q)) {
re[i]=i*choose(i-1,hit-1)*nPr(Q,i-hit)*nPr(P,hit)/nPr(N,i)
}
return(sum(re))
}
708:卵の名無しさん
18/08/30 21:08:07.53 WRrKmlGk.net
lottery_sim <- function(N,p,q,hit=1){ # N:lots, p:probability of winning
x=c(rep(1,N*p),rep(0,N*(1-p))) # q:price perlot, hit: hits required to win
win=0
try=0
while(win<hit){
n=length(x)
j=sample(n,1)
win = win + as.numeric(x[j]==1)
x=x[-j]
try=try+1
}
return(try*q)
}
lottery_sim(100,0.3,q=1,hit=3)
k=1e5
re3=replicate(k,lottery_sim(N=100,p=0.3,q=1000,hit=1))
summary(re3) ; hist(re3,col='lightblue',main='',freq=FALSE)
k=1e5
re6=replicate(k,lottery_sim(N=50,p=0.6,q=2000,hit=1))
summary(re6) ; hist(re6,col='pink', main='',freq=FALSE)
709:卵の名無しさん
18/08/31 12:51:45.85 F0nmV/5C.net
draft
√(1+x) = a + bx + cx2 + dx3 + ....
x = 0 , a=1
√(1+x) = (1 + bx + cx2 + dx3 + ....)(1 + bx + cx2 + dx3 + ....)
1 + x = 1 + 2bx + (b2+2c)x2 + (2d+2bc)x3 + ...
b=1/2
b2+2c=0, c=-1/8
2d+2bc=0,d=1/16
√(1+x) = 1+ x/2 -x2/8+ x3/16 + ....
710:卵の名無しさん
18/08/31 22:34:57.06 dQmojtaG.net
library(nleqslv)
Tetra <- function(O=c(1/2,sqrt(3)/6,sqrt(2/3)),A=c(0,0,0),B=c(1,0,0),C=c(cos(pi/3),sin(pi/3),0)){
fn <- function(x,O,A,B,C){
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO
AB=B-A
AC=C-A
c(HO%*%AB,HO%*%AC)
}
fn1 <- function(x) fn(x,O,A,B,C)
x=nleqslv::nleqslv(c(1/3,1/3),fn1)$'x'
AO=A-O
BO=B-O
CO=C-O
HO=x[1]*AO+x[2]*BO+(1-x[1]-x[2])*CO
h=sqrt(sum(HO^2))
a=sqrt(sum((B-C)^2))
b=sqrt(sum((C-A)^2))
c=sqrt(sum((A-B)^2))
s=(a+b+c)/2
base=sqrt(s*(s-a)*(s-b)*(s-c))
V=1/3*base*h
return(V)
}
Tetra()
sqrt(2)/12
711:卵の名無しさん
18/08/31 23:55:29.75 dQmojtaG.net
gcd <- function(a,b) ifelse(!a%%b,b,gcd(b,a%%b))
gcd( 349163 ,7599867)
712:卵の名無しさん
18/09/02 17:21:50.04 D0qvH3Oi.net
#include<stdio.h>
#define N 5
int b[N];
void truth_table(int k){
int i;
if(k < N){
b[k] = 1;
truth_table(k+1);
b[k] = 0;
truth_table(k+1);
}
else{
for(i = 0;i < N; i++){
printf("%s ", b[i] ? "T " : "F ");
}
printf("\n");
}
}
int main(){
truth_table(0);
return 0;
}
713:卵の名無しさん
18/09/03 19:46:52.57 oknI78Sn.net
dec2nw <- function(num, N, digit = 4){
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(r)
}
714:卵の名無しさん
18/09/05 21:52:11.29 pWY0j8lL.net
man=800
woman=200
total=man+woman
pass=100
i=0:100
rr=1.2
sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > rr*(pass-i)/woman))
715:卵の名無しさん
18/09/06 11:08:11.14 H8UAdZU9.net
(choose(5,1)*(choose(1,0)*1^5 - choose(1,1)*0^5)*(choose(4,0)*4^5 - choose(4,1)*3^5 + choose(4,2)*2^5 - choose(4,3)*1^5 + choose(4,4)*0^5)
+choose(5,2)*(choose(2,0)*2^5 - choose(2,1)*1^5 + choose(2,2)*0^5)*(choose(3,0)*3^5 - choose(3,1)*2^5 + choose(3,2)*1^5 - choose(3,3)*0^5))*2
716:卵の名無しさん
18/09/06 22:00:40.50 aj5FUH/y.net
draft
5部屋 男5 女5 定員4 空室可 混合不可
m4 m1 f(3,5) choose(5,4)*choose(5,2)*2* f(3,5)
m3 m2 f(3,5) choose(5,3)*choose(5,2)*2* f(3,5)
m3 m1 m1 f(2,5) choose(5,3)*5*4*3* f(2,5)
m2 m2 m1 f(2,5) choose(5,2)*choose(3,2)*5*4*3*f(2,5)
f(3,5)
f4 f1 f0 3!*choose(5,4)
f3 f2 f0 3!*choose(5,3)
f3 f1 f1 3!choose(5,3)
f2 f2 f1 3!*choose(5,2)*choose(3,2)
f(2,5)
f4 f1 choose(5,4)*2
f3 f2 choose(5,3)*2
717:卵の名無しさん
18/09/07 14:53:17.74 a5V/3Sw3.net
# 5部屋 男5 女5 定員4 空室可 混合不可
# [117000,] 5 5 5 5 4 3 3 3 3 2
# m4 m1 3^5-3 5P2*5
# m3 m2 3^5-3 5P2*5C2
# m3 m1 m1 2^5-2 5C3*3*5!/3!
# m2 m2 m1 2^5-2 5C3*3*5!/2!/2!
5*4*5 *240
5*4*10 *240
10*3*5*4 *30
10*3*120/2/2 *30
5*4*5*240+5*4*10*240+10*3*5*4*30+10*3*120/2/2*30
718:卵の名無しさん
18/09/08 06:48:44.13 5Jdei1ET.net
library(gtools)
n=5 # rooms
m=5 # men
w=5 # women
cap=4 # capacity
perm=permutations(n=5,r=m+w,rep=T)
nomixcap2 <- function(x,vacancy=FALSE){
men = x[1:m]
women= x[(m+1):(m+w)]
if(any(men %in% women)) return(FALSE)
if(max(tabulate(x)) > cap) return(FALSE)
if(vacancy) return(TRUE)
if(!all(1:n %in% x)) return(FALSE) # if not(all room used)
return(TRUE)
}
revac=perm[,which(apply(perm,1,function(x)nomixcap2(x,vac=T)))]
719:卵の名無しさん
18/09/08 10:20:22.08 YKQ7uLEM.net
library(gtools)
n=5 # rooms
m=5 # men
w=5 # women
perm=gtools::permutations(n,r=m+w,rep=T)
mix <- function(x){
men = x[1:m]
women= x[(m+1):(m+w)]
if(!all(men %in% women)) return(FALSE)
else all(women %in% men))
}
revac=perm[,which(apply(perm,1,mix)))]
720:卵の名無しさん
18/09/08 12:53:31.19 cGu8TtfR.net
library(gtools)
n=5 # rooms
m=5 # men
w=5 # women
perm=gtools::permutations(n,r=m+w,rep=T)
mix <- function(x){
men = x[1:m]
women= x[(m+1):(m+w)]
if(!all(men %in% women)) return(FALSE)
all(women %in% men)
}
(x=perm[sample(1:5^10,1),])
mix(x)
re=perm[which(apply(perm,1,mix)),]
head(re)
tail(re)
721:卵の名無しさん
18/09/09 07:52:51.35 C6iRFzVv.net
k=20
pair=NULL
for(i in 1:k){
a=sort(rpois(2,7))
b=sort(rpois(2,7))
c=sum(a>b)
print(c(a=a,b=b,pair=c))
pair[i]=c
}
c(pair=pair)
hist(pair)
722:卵の名無しさん
18/09/09 22:17:25.84 +VpBquyc.net
k=20
pair=NULL
for(i in 1:k){
a=sort(rpois(2,7))
b=sort(rpois(2,7))
c=sum(a>b)
print(c(a=a,b=b,pair=c))
pair[i]=c
}
hist(pair)
pair=NULL
for(i in 1:k){
a=sort(rpois(2,7))
b=sort(rpois(2,7))
c=sum(a>b)
print(c(a=a,b=b,pair=c))
pair[i]=c
}
k=10000
pair=NULL
for(i in 1:k){
a=sort(jitter(rpois(2,7)))
b=sort(jitter(rpois(2,7)))
c=sum(a>b)
pair[i]=c
}
hist(pair,col='skyblue')
723:卵の名無しさん
18/09/09 22:36:16.30 +VpBquyc.net
# reverse' [] = []
# reverse' (x:xs) = reverse'(xs) ++ [x]
#
# main = do
# print $ reverse' [1..10]
reverse <- function(x){
if(!length(x)) return(NULL)
c(Recall(x[-1]),x[1])
}
cat(reverse(LETTERS[1:26]))
724:卵の名無しさん
18/09/10 14:53:39.38 40BngIEI.net
draft
rot13ch <- function(ch){
if('A' <= ch & ch <= 'M' ) LETTERS[which(x==LETTERS)+13]
if('a' <= ch & ch <= 'm' ) LETTERS[which(x==letters)+13]
if('N' <= ch & ch <= 'Z' ) LETTERS[which(x==LETTERS)-13]
if('n' <= ch & ch <= 'z' ) LETTERS[which(x==letters)-13]
else ch
}
rot13 <- function(x){
x=strsplit(x,NULL)
re=NULL
while(!length(x)){
re=append(re,rot13ch(x[1]))
x=x[-1]}
return(re)
}
(hello13 = rot13 ("Uraguchi!"))
print(rot13( hello13))
725:卵の名無しさん
18/09/10 15:06:30.76 40BngIEI.net
draft2
rot13ch<- function(x){
if('A' <= x & x <= 'M' ) re=LETTERS[which(x==LETTERS)+13]
if('a' <= x & x <= 'm' ) re=letters[which(x==letters)+13]
if('N' <= x & x <= 'Z' ) re=LETTERS[which(x==LETTERS)-13]
if('n' <= x & x <= 'z' ) re=letters[which(x==letters)-13]
else re=x
return(re)
}
rot13ch('j')
rot13 <- function(x){
x=strsplit(x,NULL)
re=NULL
while(!length(x)){
re=append(re,rot13x(x[1]))
x=x[-1]}
return(re)
}
726:卵の名無しさん
18/09/10 16:11:10.56 40BngIEI.net
draft 3
rot13ch<- function(x){
if('A' <= x & x <= 'M' ) re=LETTERS[which(x==LETTERS)+13]
else
if('a' <= x & x <= 'm' ) re=letters[which(x==letters)+13]
else
if('N' <= x & x <= 'Z' ) re=LETTERS[which(x==LETTERS)-13]
else
if('n' <= x & x <= 'z' ) re=letters[which(x==letters)-13]
else re=x
return(re)
}
rot13 <- function(str){
x=unlist(strsplit(str,NULL))
re=NULL
while(!length(x)){
re=append(re,rot13ch(x[1]))
x=x[-1]}
return(re)
}
str="Uraguchi!"
rot13(str)
727:卵の名無しさん
18/09/10 18:57:39.99 V03Ln2wY.net
Gacha <- function(p){ # p: probability of each Gacha item
p=p/sum(p)
sum.rev <- function(x){ # i,j,k -> 1/(p[i]+p[j]+p[k])
n=length(x)
s=numeric(n)
for(i in 1:n) s[i]=p[x[i]]
1/sum(s)
}
n=length(p)
re=numeric(n)
for(i in 1:n) re[i]=(-1)^(i-1)*sum(apply(combn(n,i),2,sum.rev))
sum(re)
}
Gacha(1:4/10)
Gacha(c(1/10,rep(9/50,5)))
#
sim <- function(p){
p=p/sum(p)
n=length(p)
y=NULL
while(!all(1:n %in% y)){
y=append(y,sample(1:n,1,prob=p))
}
return(length(y))
}
mean(replicate(1e4,sim(1:4)))
Gacha(1:4)
mean(replicate(1e5,sim(c(9,9,9,9,9,5))))
Gacha(c(9,9,9,9,9,5))
728:卵の名無しさん
18/09/10 20:37:33.27 YF8r2q2P.net
１万回のシミュレーション
> mean(replicate(1e4,sim(1:10)))
[1] 68.7788
理論値
> Gacha(1:10)
[1] 68.98458
729:卵の名無しさん
18/09/11 13:14:50.53 DPX9K4Dn.net
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]]
re1
for(i in 2:(n-1)){
re1=c(re1,ifelse(i%%2,' + ',' - '),'(',re[[i]],')')
}
output=c(paste(re1,collapse=""),ifelse(n%%2,'+','-'), re[[n]])
cat(output,'\n')
if(write) write(output,'output.txt')
invisible(output)
}
Gacha.fm(c(1/10,2/10,3/10,4/10))
Gacha.fm(c(9/50,9/50,9/50,9/50,9/50,5/10))
730:卵の名無しさん
18/09/11 13:15:42.26 DPX9K4Dn.net
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]]
re1
for(i in 2:(n-1)){
re1=c(re1,ifelse(i%%2,' + ',' - '),'(',re[[i]],')')
}
output=c(paste(re1,collapse=""),ifelse(n%%2,'+','-'), re[[n]])
cat(output,'\n')
if(write) write(output,'output.txt')
invisible(output)
}
Gacha.fm(c(1/10,2/10,3/10,4/10))
Gacha.fm(c(9/50,9/50,9/50,9/50,9/50,5/10))
731:卵の名無しさん
18/09/12 17:03:31.45 +rWfu75u.net
draft
insert <- function(x,yys){
if(!length(yys)) return x
else {
y=yy[1]
ys=yys[-1]
if(x<y) return c(x,y,ys)
else c(y, Recall(x,ys))
}
}
isort <- function(xxs){
if(!length(xxs)) return NULL
else{
x=xxs[1]
xs=xxs[-1]
insert(x,Recall(xs))
}
isort(c(4, 6, 9, 8, 3, 5, 1, 7, 2))
732:卵の名無しさん
18/09/12 17:41:27.26 MBeBhzv/.net
insert <- function(x,y){
if(!length(y)) return(x)
if(x<y[1]) return(c(x,y))
return(c(y[1], Recall(x,y[-1])))
}
isort <- function(x){
if(!length(x)) return(NULL)
insert(x[1],Recall(x[-1]))
}
isort(c(4, 6, 9, 8, 3, 5, 1, 7, 2))
733:卵の名無しさん
18/09/12 20:36:00.82 +rWfu75u.net
bsort <- function(x){
if(!length(x)) return(NULL)
if(length(x)==1) return(x)
else{
y=bsort(x[-1])
ifelse(x[1]<y[1],c(x[1],y),c(y[1],bsort(c(x[1],y[-1])
734:} } bsort(c(3,1,4,1,5,9,2,6,5)
735:卵の名無しさん
18/09/13 07:37:31.93 ha3KPNlb.net
sim2 <- function(p){
n=length(p) # number of items
if(sum(p)>=1){ # no blank and/or rate of probabilities
prob=p/sum(p) # scaling for sum(prob)=1
lot=1:n # no blank lot
}else{
prob=c(p,1-sum(p)) # blank with probability of 1-sum(p)
lot=1:(n+1) # lot[n+1] blank lot
}
y=NULL
count=0
while(length(y)<n){
z=sample(lot,1,prob=prob)
count=count+1
if(!any(z==y)) y=append(y,z) # append new item only
}
return(count)
}
736:卵の名無しさん
18/09/13 08:15:32.46 ha3KPNlb.net
p=c(9,9,9,9,9,5)/50
n=length(p)
seg=cumsum(p)
p
seg
x=runif(1)
x
x > seg
k=numeric(n)
k | (x > seg)
737:卵の名無しさん
18/09/13 08:35:57.43 ha3KPNlb.net
p=c(9,9,9,9,9,5)/50
n=length(p)
seg=cumsum(p)
count=0
k=numeric(n)
while(k<2^n-1){
x=runif(1)
k | (x > seg)
k=k | (x > seg)
count=count+1
}
count
738:卵の名無しさん
18/09/13 11:52:48.72 vlGxfg9U.net
sim1 <- function(p){
n=length(p)
sep=cumsum(p)
y=NULL
count=0
while(length(y) < n){
z=sample(1:n,1,prob=p)
if(!any(z==y)) y=append(y,z) # append new item only
count=count+1
}
return(count)
}
sim3 <- function(p){
n=length(p)
sep=cumsum(p)
y=NULL
count=0
while(length(y) < n){
z=sum(runif(1) < sep)
if(!any(z==y)) y=append(y,z) # append new item only
count=count+1
}
return(count)
}
p=c(9,9,9,9,9,5)/50
system.time(mean(replicate(1e4, sim1(p))))
system.time(mean(replicate(1e4, sim3(p))))
739:卵の名無しさん
18/09/13 19:14:06.79 ha3KPNlb.net
qsort [] = []
qsort (n:xs) = qsort lt ++ [n] ++ qsort gteq
where
lt = [x | x <- xs, x < n]
gteq = [x | x <- xs, x >= n]
main = do print $ qsort [4, 6, 9, 8, 3, 5, 1, 7, 2]
740:卵の名無しさん
18/09/14 07:47:53.11 IOEu7JkF.net
foo <- function(a,b,c){
gr=expand.grid(1:a,1:b)
c2=(1:c)^2
f = function(x,y) (x<y) & (x^2+y^2) %in% c2
i=which(mapply(f,gr[,1],gr[,2]))
ab=as.matrix(gr[i,])
cbind(ab,sqrt(ab[,1]^2+ab[,2]^2))
}
foo(20,20,20)
a=20;b=20;c=20
gr=expand.grid(1:a,1:b,1:c)
f = function(a,b,c) a<b & a^2+b^2==c^2
i=which(mapply(f,gr[,1],gr[,2],gr[,3]))
gr[i,]
741:卵の名無しさん
18/09/14 08:01:47.72 IOEu7JkF.net
foo <- function(a,b,c){
gr=expand.grid(1:a,1:b)
c2=(1:c)^2
f = function(x,y) (x<y) & (x^2+y^2) %in% c2
i=which(mapply(f,gr[,1],gr[,2]))
ab=as.matrix(gr[i,])
cbind(ab,sqrt(ab[,1]^2+ab[,2]^2))
}
foo(20,20,20)
foo(100,100,100)
a=20;b=20;c=20
fooo <- function(a,b,c){
gr=expand.grid(1:a,1:b,1:c)
f = function(a,b,c) a<b & a^2+b^2==c^2
i=which(mapply(f,gr[,1],gr[,2],gr[,3]))
gr[i,]
}
fooo(20,20,20)
fooo(100,100,100)
742:卵の名無しさん
18/09/14 09:51:33.26 pDMF5zq8.net
#include<stdio.h>
#include<stdlib.h>
int pit(int A,int B, int C){
int a,b,c; int pit=0;
for(a=1;a<=A;a++){
for(b=a;b<=B;b++){
for(c=b;c<=C;c++){
if((a<b) && (a*a+b*b==c*c)) {
++pit;
printf("%d : %d * %d + %d * %d = %d * %d\n",pit,a,a,b,b,c,c);
}
}
}
}
return 0;
}
int main( int argc, char *argv[] ){
int a,b,c;
a = atoi(argv[1]);
b = atoi(argv[2]);
c = atoi(argv[3]);
pit(a,b,c);
return 0;
}
743:卵の名無しさん
18/09/15 09:37:12.07 LVUNnMZD.net
とある会社の社長は毎日午後5時に会社を出て自宅からの迎えのクルマに乗って帰る。
ある日、午後4時に退社した。
天気が良かったので、迎えのクルマに出会うまで散歩した。
出会ったところで、クルマはUターンして自宅に戻った。
するといつもより10分早く帰宅した。
何時何分にクルマに出会ったか？
URLﾘﾝｸ(cybozushiki.cybozu.co.jp)
尚、迎えの車は５時に会社に到着するように自宅を出発し行きも帰りも等速走行を仮定する。
744:卵の名無しさん
18/09/15 18:48:34.28 LVUNnMZD.net
診療所から病院に患者を救急搬送する。
病院から救急車が診療所に向かっており１０時到着予定と連絡が入った。
患者が９時に急変したため診療所の普通車で病院に向かって救急車と出会ったら救急車に患者を移して
745:搬送し病院到着を早めることになった。当然、救急車の方が速く走れる。 ９時５０分に救急車に乗り移ることができた。 病院到着は予定より何分早まるか述べよ。 乗り換えに要する時間は０とする。
746:卵の名無しさん
18/09/15 22:29:13.55 LVUNnMZD.net
>>490
一般化してみた。
歩行時間:t
歩行距離:l
迎えの車速:rv
同乗の車速:v
出発時刻差:s
到着時刻差:X
通常走行時間:d/rv+d/v
早退時走行時間:(d-l)/rv+(d-l)/v
d/rv+d/v - ((d-l)/rv+(d-l)/v)= X
l (r + 1)/rv=X
l/v=X*r/(r+1)
s0+s+d/v - {s0+t+(d-l)/v}= X
s-t+l/v=X
X=(r + 1) (s - t)
t = s - X / (r+1)
747:卵の名無しさん
18/09/16 18:25:50.55 MIF83oNx.net
# include<stdio.h>
# include <stdlib.h>
# define N 10
int a[N + 1], key;
int search(int a[],int key,int i);
int main(int argc, char *argv[])
{
int a[]= {0,1,2,3,4,5,6,7,8,9};
key = atoi(argv[1]);
a[N]=key;
if (search(a, key ,0) < N) {
printf("Found!\n");
}
else{
printf("Not Found!\n");
}
}
int search(int a[], int key, int i)
{
if(a[i]==key){
return i;
}
return search(a,key,++i);
}
748:卵の名無しさん
18/09/16 19:43:53.55 MIF83oNx.net
# 歩行時間:t
# 歩行速度:w
# 迎えの車速:v0
# 同乗の車速:v1
# 出発時刻差:s
# 到着時刻差:X
# l=wt
# 通常走行時間:d/v0+d/v1
# 早退時走行時間:(d-wt)/v0+(d-wt)/v1
# d/v0+d/v1 - ((d-wt)/v0+(d-wt)/v1)= X
# wt(1/v0+1/v1) = X
#
# s0+s+d/v1 - {s0+t+(d-wt)/v1}= X
# s - t(1-w/v1) = X
#
s = t + wt/v0
t = sv0/(v0+w)
X = tw(1/v0+1/v1) = sv0/(v0+w)*w(1/v0+1/v1)
749:卵の名無しさん
18/09/16 21:14:22.32 MIF83oNx.net
# The ambulance arrives X hours ealier at hospital,
# when the patient leave clinic s hours earlier than planned
# and encounter the ambulance t hours later,
# ambulance runs with the velocity of v0 without patient, and v1 with patient
# clinic car runs with the velocity of w
Earlier <- # s hour earlier departure, X hour earlier arrival
function(s=NULL,X=NULL,t=NULL,v0=60,v1=45,w=30){
if(is.null(s)) re=c(s=X/(v0/(v0+w)*w*(1/v0+1/v1)),t=X/w/(1/v0+1/v1))
if(is.null(X)) re=c(X=s*v0/(v0+w)*w*(1/v0+1/v1),t=s*v0/(v0+w))
if(!is.null(t)) re=c(s=t + t*w/v0, X = t*w*(1/v0+1/v1))
return(re)
}
Earlier(X=10/60)*60
Earlier(s=30/60)*60
Earlier(t=30/60)*60
750:卵の名無しさん
18/09/17 08:24:17.97 Oiy+BYJP.net
こういう計算ができるとdoor to balloon timeが短縮できるから臨床医に必要な能力だな。
診療所から病院に患者を救急搬送する。
病院から医師搭乗の救急車が診療所に向かっており１０時到着予定と連絡が入った。
患者の病態が悪化したら、診療所の普通車で病院に向かい救急車と出会ったら
救急車に患者を移して搬送し病院到着を急ぐという計画を立てた。
普通車から救急車への患者の乗り換えで１０分余分に時間がかかる。
道路事情から病院から診療所への道は
平均時速６０ｋｍで、逆方向は平均時速４５ｋｍで定速走行する。診療所の普通車は信号待ちもあり平均時速３０kmで定速走行する。
何時以降の病態悪化は診療所の車を使わずに救急車の到着を待つ方が病院に早く着くか？
751:卵の名無しさん
18/09/17 08:25:24.95 Oiy+BYJP.net
こういう計算ができるとdoor to balloon timeが短縮できるから臨床医に必要な能力だな。
診療所から病院に患者を救急搬送する。
病院から医師搭乗の救急車が診療所に向かっており１０時到着予定と連絡が入った。
患者の病態が悪化したら、診療所の普通車で病院に向かい救急車と出会ったら
救急車に患者を移して搬送し病院到着を急ぐという計画を立てた。
普通車から救急車への患者の乗り換えで１０分余分に時間がかかる。
道路事情から救急車は病院から診療所への道は
平均時速６０ｋｍで、逆方向は平均時速４５ｋｍで定速走行する。診療所の普通車は信号待ちもあり平均時速３０kmで定速走行する。
何時以降の病態悪化は診療所の車を使わずに救急車の到着を待つ方が病院に早く着くか？
752:卵の名無しさん
18/09/17 20:23:03.24 cb+FssaI.net
診療所から病院に患者を救急搬送する。
病院から救急車が診療所に向かっており１０時到着予定と連絡が入った。
患者が９時に急変したため診療所の普通車で病院に向かって救急車と出会ったら救急車に患者を移して搬送し病院到着を早めることになった。救急車の方が速く走れる。
９時５０分に救急車に乗り移ることができた。
病院到着は予定より何分早まるか述べよ。
車は定速走行とし、乗り換えに要する時間は考慮しない。
753:卵の名無しさん
2018/09/1
754:7(月) 22:45:36.98 ID:cb+FssaI.net
755:卵の名無しさん
18/09/17 22:46:00.21 cb+FssaI.net
# while loop version for dec2n
dec2nw <- function(num, N, digit = 4){
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(r)
}
756:卵の名無しさん
18/09/18 01:40:12.09 jA0T/PMd.net
URLﾘﾝｸ(www.tutorialspoint.com)
757:卵の名無しさん
18/09/18 06:21:30.53 jA0T/PMd.net
>>701
URLﾘﾝｸ(tpcg.io)
758:卵の名無しさん
18/09/18 07:43:40.10 jA0T/PMd.net
#include<stdlib.h>
#define N 5
int b[N];
//char * u[] = {"0","1","2","3","4"};
char * u[] = {"志学","而立","不惑","知命","耳順","従心"};
void fill(int k){
int i,j;
if(k < N){
for(i=0;i<N;i++){
for(j=0; j < k && b[j]!=i; j++);
;
if(j==k){
b[k] = i;
fill(k+1);
}
}
}
else{
for(i = 0;i < N; i++){
printf("%s ", u[b[i]] );
}
printf("\n");
}
}
int main(){
printf("all around\n");
fill(0);
return 0;
}
759:卵の名無しさん
18/09/19 22:38:29.64 Eiqd7mE7.net
draft
is.sorted <- function(x){
if(length(x)==1) return(TRUE){
else{
if(x[1]<=x[2]) is.sort(x[-1]){
else return(FALSE)
}
}
}
760:卵の名無しさん
18/09/19 22:44:21.96 QK/Jo93v.net
is.sorted <- function(x){
if(length(x)==1){ return(TRUE)
}else{
if(x[1]<=x[2]) is.sorted(x[-1])
else return(FALSE)
}
}
761:卵の名無しさん
18/09/20 20:45:13.64 iA1AjZuZ.net
ニューロン治療できる医療大麻オイルの紹介
URLﾘﾝｸ(plaza.rakuten.co.jp)
ニューロン=人工知能のモデルとなっている神経細胞のやりとりするところ
この細胞の伝達する場所に大麻受容体があります この受容体を通って治療効果を得る大麻を医療大麻と呼びます
人間に大麻受容体があったなんて不思議ですよね
762:卵の名無しさん
18/09/22 08:25:11.10 5EtHM4o7.net
>>711 は国試に２３回も落ちてて
くるくるぱーの裏口バカに
なっちゃってるのらぁあぁぁ
Ｆラン事務員の濃ゆぅぅい生ガキ汁
ド底辺の臭いが落ちないよぉ
んほおぉぉぉおぉぉ
763:卵の名無しさん
18/09/22 08:33:45.20 zrIlwkkB.net
>>713
>711って再帰呼び出ししているけど何をするスクリプトかわかる？
764:卵の名無しさん
18/09/22 13:39:33.84 NI3pjyOo.net
#include<stdio.h>
#include<stdlib.h>
#define N 8
#define W 2*N-1
int b[N],row[N],up[W],down[W],count=0;
void fill(int k){
int i,j;
if(k < N){
for(i=0;i<N;i++){
if(!row[i] && !up[i+k] && !down[i-k+N-1]){
b[k] = i;
row[i]=up[i+k]=down[i-k+N-1]=1;
fill(k+1);
row[i]=up[i+k]=down[i-k+N-1]=0;
}
}
}
else{
printf("# %2d\n",++count);
for(i =0; i < N; i++){
for(j = 0; j < N; j++){
if(b[i] == j) printf(" Ｑ ");
else printf(" ・ ");
}
printf("\n");
}
}
}
int main(){
fill(0);
return 0; }
765:卵の名無しさん
18/09/22 21:34:00.39 zrIlwkkB.net
#
766:include <stdio.h> int arr[] ={1,2,3,4,5}; int total = (sizeof(arr) / sizeof(arr[0])); void heapify(int arr[], int i) { int lft = i * 2; int rgt = lft + 1; int grt = i; if (lft <= total && arr[lft] > arr[grt]) grt = lft; if (rgt <= total && arr[rgt] > arr[grt]) grt = rgt; if (grt != i) { int tmp=arr[grt]; arr[grt]=arr[i]; arr[i]=tmp; heapify(arr, grt); } } int main(){ heapify(arr,0); for (int i = 0; i < total; i++) printf("%d ", arr[i]); printf("\n"); }
767:卵の名無しさん
18/09/23 18:58:09.75 edPyfcCi.net
library(gtools)
k=3
perm=unique(permutations(3*k,3*k,rep(1:3,k),set=F,rep=F))
is3 <- function(x){
n=length(x)
y=c(x,x[1])
re=NULL
for(i in 1:(n-1)) re[i]=all((1:3) %in% y[i:(i+2)])
any(re)
}
i=(!apply(perm, 1, is3))
perm1=perm[i,]
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))
}
768:卵の名無しさん
18/09/23 18:58:40.06 edPyfcCi.net
same <- function(x,base){
f=function(a,b=base){ # is equal to base
all(a==b)
}
mat=rbind(rn(x),rn(rev(x))) # rotation and/or symmetric conversion
any(apply(mat,1,f))
}
mat=perm1
tail(perm1)
exvar <- function(row,mat=perm1){
var=rn(mat[row,])
n.var=nrow(var)
idx=NULL
for(i in 1:n.var){
idx=c(idx,which(apply(mat,1,function(x) same(x,var[i,]))))
}
min(idx)
}
index=sapply(1:nrow(perm1),exvar) # several minutes for k=3
re=perm1[unique(index),]
re
print(matrix(c("赤","白","青")[re],ncol=3*k),quote=F)
769:卵の名無しさん
18/09/24 16:11:21.71 cwdsebwG.net
same <- function(x,base){
f=function(a,b=base){ # is equal to base
all(a==b)
}
mat=rbind(rn(x),rn(rev(x))) # rotation and/or symmetric conversion
any(apply(mat,1,f))
}
mat=perm1
core=NULL # 基本解を集める行列
while(nrow(mat)){
var=rn(mat[1,]) # matの先頭行のバリアント
n.var=nrow(var)
idx=NULL # 先頭行とそのバリアントのindex
for(i in 1:n.var){ # 先頭行およびバリアント行のindexを返す
idx=c(idx,which(apply(mat,1,function(x) same(x,var[i,]))))
}
i=sort(unique(idx))
core=rbind(core,mat[i[1],])
mat=mat[-i,] # 一度、検討した行と同一相当行を除いたmatで再評価に回す
}
core
print(matrix(c("赤","白","青")[core],ncol=3*k),quote=F)
770:卵の名無しさん
18/09/24 20:26:43.18 cwdsebwG.net
library(compiler)
cuniq=cmpfun(unique)
jperm <- function(n,r,v) {
if (r == 1)
matrix(v, n, 1)
else if (n == 1)
matrix(v, 1, r)
else {
X <- NULL
for (i in 1:n) X <- rbind(X, cbind(v[i], jperm(n - 1, r - 1, v[-i])))
cuniq(X)
}
}
cjperm=cmpfun(jperm)
v=rep(1:3,4)
n=r=length(v)
perm4=cjperm(n,r,v)
771:卵の名無しさん
18/09/24 20:44:38.00 cwdsebwG.net
#include<stdio.h>
#define N 6
int b[N];
char * u[] = {"P","E","P","P","E","R"};
void p_table(int k){
int i;
if(k < N){
for(i=0;i<N;i++){
b[k] = i;
p_table(k+1);
}
}
else{
for(i = 0;i < N; i++){
printf("%s ", u[b[i]] );
}
printf("\n");
}
}
int main(){
p_table(0);
return 0;
}
772:卵の名無しさん
18/09/24 21:05:06.98 cwdsebwG.net
// URLﾘﾝｸ(www.geeksforgeeks.org)
// Program to print all permutations of a string in sorted order.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
/* Following function is needed for library function qsort(). */
int compare(const void *a, const void * b)
{
return ( *(char *)a - *(char *)b );
}
// A utility function two swap two characters a and b
void swap(char* a, char* b)
{
char t = *a;
*a = *b;
*b = t;
}
773:卵の名無しさん
18/09/24 21:05:27.97 cwdsebwG.net
// This function finds the index of the smallest character
774: // which is greater than 'first' and is present in str[l..h] int findCeil(char str[], char first, int l, int h) { // initialize index of ceiling element int ceilIndex = l; // Now iterate through rest of the elements and find // the smallest character greater than 'first' for (int i = l+1; i <= h; i++) if (str[i] > first && str[i] < str[ceilIndex]) ceilIndex = i; return ceilIndex; }
775:卵の名無しさん
18/09/24 21:05:44.12 cwdsebwG.net
// Print all permutations of str in sorted order
void sortedPermutations(char str[])
{
// Get size of string
int size = strlen(str);
// Sort the string in increasing order
qsort(str, size, sizeof( str[0] ), compare);
// Print permutations one by one
bool isFinished = false;
while (!isFinished)
{
// print this permutation
static int x = 1;
printf("%d %s \n", x++, str);
// Find the rightmost character which is smaller than its next
// character. Let us call it 'first char'
int i;
for (i = size - 2; i >= 0; --i)
if (str[i] < str[i+1])
break;
776:卵の名無しさん
18/09/24 21:06:30.92 cwdsebwG.net
// If there is no such chracter, all are sorted in decreasing order,
// means we just printed the last permutation and we are done.
if (i == -1)
isFinished = true;
else
{
// Find the ceil of 'first char' in right of first character.
// Ceil of a character is the smallest character greater than it
int ceilIndex = findCeil(str, str[i], i + 1, size - 1);
// Swap first and second characters
swap(&str[i], &str[ceilIndex]);
// Sort the string on right of 'first char'
qsort(str + i + 1, size - i - 1, sizeof(str[0]), compare);
}
}
}
// Driver program to test above function
int main()
{
char str[] = "";
sortedPermutations( str );
return 0;
}
777:卵の名無しさん
18/09/26 11:07:08.91 Rj8Wp2L9.net
b=0
w=1
bwb=function(x){
n=length(x)
if(n<3) return(F)
else x[n-2]==b & x[n-1]==w & x[n]==b
}
f = function(){
x=NULL
while(!bwb(x)){
x=c(x,sample(c(b,w),1,prob=c(0.2,0.8)))
}
c(W=sum(x==w),B=sum(x==b))
}
re=replicate(1e4,f())
mean(re[1,])
mean(re[2,])
778:卵の名無しさん
18/09/26 22:08:29.29 +zbP8/80.net
library(Rmpfr)
f <- function(NN,prec=1000){
for(k in 1:NN){
for(m in k:NN){
for(n in m:NN){
K=mpfr(factorial(k),prec)
M=mpfr(factorial(m),prec)
N=mpfr(factorial(n),prec)
if(K+M==N) cat(k,'! + ',m,'! = ',n,'!\n')
}
}
}
}
try(f(100))
779:卵の名無しさん
18/10/02 21:03:43.40 DMIGI67G.net
print $ foldr (-) 0 [1,2,3]
print $ foldr (-) 0 ( 1:(2:(3:[])) )
print $ 1-(2-(3-0))
780:卵の名無しさん
18/10/02 22:17:49.03 DMIGI67G.net
elem_r y ys = foldr (\x acc -> if x==y then True else acc) False ys
elem_r 39 [3.6..]
elem_l y ys = foldl(\acc x -> if x==y then True else acc) False ys
{- ERROR!!
elem_l 39 [3,6..]
781:卵の名無しさん
18/10/03 13:45:16.24 OGYXL17c.net
map' f xs = foldr (\x acc -> f x:acc) [] xs
elem_r y ys = foldr (\x acc -> if x==y then True else acc) False ys
elem_l y ys = foldl(\acc x -> if x==y then True else acc) False ys
main = do
print $ zipWith (\x y -> (x*30+3)/y) [5,4,3,2,1] [1..5]
print $ map (\(a, b) -> a+b) [(1,2),(3,5),(6,3),(2,6),(2,5)]
print $ zipWith (flip (++)) ["love you","is an angel"] ["I ","Mei "]
print $ map (flip subtract 20) [1..5]
print $ foldl (\a b -> a+b) 0 [3,5,2,1]
print $ foldl (-) 0 [1,2,3]
print $ ((0-1)-2)-3
print $ foldr max 0 [1,2,3]
print $ 1 `max` (2 `max`(3 `max` 0))
print $ foldr (-) 0 [1,2,3]
print $ foldr (-) 0 ( 1:(2:(3:[])) )
print $ 1-(2-(3-0))
print $ foldr (\a b -> a-b) 0 ( 1:(2:(3:[])) )
782:卵の名無しさん
18/10/08 17:55:51.64 pH+BqnrA.net
import Data.List
divisor n = find (\m -> n `mod` m ==0 )[2..floor.sqrt.fromIntegral $ n]
main = do
print $ divisor $ (10^19-1) `div` 9
783:卵の名無しさん
18/10/08 18:06:39.22 pH+BqnrA.net
import Data.List
divisor n = find (\m -> n `mod` m ==0 )[2..floor.sqrt.fromIntegral $ n]
divisors n = filter (\m -> n `mod` m ==0 )[2..n-1]
main = do
print $ divisor $ (10^19-1) `div` 9
784:卵の名無しさん
18/10/08 20:25:22.38 pH+BqnrA.net
１００万個の○が円形に並んでいます。
図のように、まず１つの○に色をぬり、
次にその●から時計回りに１０８個進んで止まり、そこにある○をぬります。
さらにその●から
時計回りに１０８個進んで止まり、
そこにある○をぬり、以下同じ作業を くり返していきます。
すでに色がぬられた●に止まったときに終了とするとき、
何個の○をぬることができますか？
foo jump spot = length $ takeWhile (\x -> x*jump `mod` spot /=0) [1..]
1 + foo 108 1000000
785:卵の名無しさん
18/10/08 20:30:18.22 pH+BqnrA.net
ｍ個の○が円形に並んでいます。
図のように、まず１つの○に色をぬり、
次にその●から時計回りにｎ個進んで止まり、そこにある○をぬります。
さらにその●から
時計回りにｎ個進んで止まり、
そこにある○をぬり、以下同じ作業を くり返していきます。
すでに色がぬられた●に止まったときに終了とするとき、
何個の○をぬることができますか？
ｍ / gcd(m.n)
m,nが互いに素なら全部塗れる。
786:卵の名無しさん
18/10/08 20:48:49.62 pH+BqnrA.net
gcd a b
| a `mod` b == 0 = b
| otherwise = gcd b $ a `mod` b
787:卵の名無しさん
18/10/10 08:15:31.85 VVThWPPb.net
事務員 とは？
同一者の別名として国試浪人、統計野郎、自称医科歯科、がある
自称医科歯科卒、実際は九州の駅弁国立出身で、卒業後は実家の東京に帰り国試浪人となり23年間経過、家庭教師バイトなどを経て現在は事務員、とされている
本人は、勤務医でたまに開業医の手伝いと内視鏡バイト、専門医なし、と主張しているが、連日連夜の異常な書き込み数の多さからは、勤務医とは考え難い
彼の職業がたとえ何であったとしても、人としてド底辺なことは間違いがない
自ら事務員であると告白したこともある
ｽﾚﾘﾝｸ(hosp板:108番)-109
自作自演も頻繁に行なっている、内視鏡バイトの後はステーキハウスに行く妄想話がよく出てくる、実際には食べたこともないんだろうな
病名を挙げて架空の診療報告を行うこともあるが、今どきヒヨッコ研修医でもそんなことやらねーぞW
事務仕事の際に手にした他人の給与明細や同窓会費振込票を盗撮し、自分のものとしてアップロードしたこともある、犯罪なんじゃねーの？
統計絶対主義でマウントを取ろうとするが、大学卒業後の人生の大部分においてそれが何の役に立つのかはなはだ疑問である
統計に基づいた重病の治療成績を患者に説明しても、患者の希望する治療は患者の事情によって異なることは一定年数以上の臨床経験があるものなら誰でも知っていることだ
それが理解できない彼は、彼自身が素性を明らかにしているようなものである
コンプ薬屋同様に、私立医への執拗な誹謗中傷を数年余にわたって続けている
裏口批判をしているが、実際にそれが可能な資産家の子弟がいるであろう右派スレには近寄らない、一度書き込�
788:ｾ形跡があったが指摘したら別人だと言って逃げた 弱いものが多そうなスレしか狙わないのもコンプ薬屋と共通する異常者の特徴である 劣等コンプレックスの塊なんだろうな 自分の主張に対して都合の悪い突込みが入ると、相手を私立医と決めつけて、さてはシリツだな、裏口バカ、というセリフでワンパターンに返すことが多い 国試問題を挙げて簡単だとほざいてますが、本番で通らなければお話になりません いくらネットでわめいていても、医師国試の本番で通らなければお話になりません 大事なことなので２回言いました
789:卵の名無しさん
18/10/12 06:12:09.09 co72L/r8.net
知り合いから教えてもらった自宅で稼げる方法
興味がある人はどうぞ
みんながんばろうねぇ『羽山のサユレイザ』で
K5K
790:卵の名無しさん
18/10/12 20:06:27.81 s0yH21ZM.net
ある医大で合格率の男女比が1.2で男子有意という結果だったという。
定員100で男子800人女子200人が受験して合格率の男女比が
1.2であったときに統計的には有意差があると言えるか？
791:卵の名無しさん
18/10/13 08:03:38.90 cNomfA2G.net
man=900
woman=100
pass=100
total=man+woman
p=pass/total
i=0:100
#
sum((choose(woman,i)*choose(man,pass-i)/choose(total,pass))*(i/woman > (pass-i)/man))
sum((choose(woman,i)*choose(man,pass-i)/choose(total,pass))*(i/woman == (pass-i)/man))
sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > (pass-i)/woman))
sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > 1.2*(pass-i)/woman))
792:卵の名無しさん
18/10/13 08:04:10.37 cNomfA2G.net
> sum((choose(woman,i)*choose(man,pass-i)/choose(total,pass))*(i/woman > (pass-i)/man))
[1] 0.4160339
> sum((choose(woman,i)*choose(man,pass-i)/choose(total,pass))*(i/woman == (pass-i)/man))
[1] 0.1389853
> sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > (pass-i)/woman))
[1] 0.4449808
> sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > 1.2*(pass-i)/woman))
[1] 0.3090334
793:卵の名無しさん
18/10/13 08:28:07.92 cNomfA2G.net
ある医大で合格率の男女比が1.2で男子有意という結果だったという。
定員100で男子800人女子200人が受験して合格率の男女比が
1.2であったときに統計的には有意差があると言えるか？
f <- function(ratio=1.2,man=800,woman=200,pass=100){
total=man+woman
i=0:pass
sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > ratio*(pass-i)/woman))
}
f()
794:卵の名無しさん
18/10/13 10:57:49.86 WXqDSsDl.net
draft
{-ある医大で合格率の男女比が1.2で男子有意という結果だったという。
定員100で男子800人女子200人が受験して合格率の男女比が
1.2であったときに統計的には有意差があると言えるか？
-}
choose n r = [1..n] `div` [1..n-r] `div` [1..r]
ratio=1.2
pass=100
man=800
woman=200
total=man+woman
ii=[0..pass]
male_fem i =
| choose(man,i)*choose(woman,pass-i)/choose(total,pass) <ratio 0
| otherwise choose(man,i)*choose(woman,pass-i)/choose(total,pass
sum $ map male_fem ii
795:卵の名無しさん
18/10/13 12:59:15.69 PvPTvQgs.net
>>736
ド底辺シリツ医大卒って医療従事者にマウントするのは不可能だから事務員や国試浪人認定するしかないんだなぁ。
こういうのにサクッと答えてド底辺シリツ医大卒でも高卒の基礎学力くらいあるのを示せばいいのに。
ある医大で合格率の男女比が1.2で男子優位という結果だったという。
定員100で男子800人女子200人が受験して合格率の男女比が
1.2であったときに統計的には有意差があると言えるか？
796:卵の名無しさん
18/10/13 15:21:07.70 PvPTvQgs.net
>>743
Rで数えると
f <- function(ratio=1.2,man=800,woman=200,pass=100){
total=man+woman
i=0:pass
sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > ratio*(pass-i)/woman))
}
f()
Haskellで数えると
ratio=1.2
pass=100
man=800
woman=200
total=man+woman
ii=[0..pass]
male_fem i= if r < ratio then 0 else p
where a = (choose man i)*(choose woman (pass-i))
b = choose total pass
m=fromIntegral i/ fromIntegral man
w=fromIntegral (pass-i)/fromIntegral woman
r = m/w
p = (fromIntegral a) / (fromIntegral b)
main = do
print $ sum ( map male_fem ii)
797:卵の名無しさん
18/10/14 06:36:53.74 +n7lCedi.net
f <- function(ratio=1.2,wm=0.2,total=1000,pass=100){
woman=total*wm
man=total-woman
i=0:pass
sum((choose(man,i)*choose(woman,pass-i)/choose(total,pass))*(i/man > ratio*(pass-i)/woman))
}
f()
r=seq(1,3,length=1000)
p=sapply(r,f)
plot(r,p,pch=19,type='n',bty='l',xlab='M/F ratio',ylab='probability')
abline(h=0.01,lty=3)
g <- function(wm) lines(r,sapply(r,function(r)f(r,wm)),col=10*wm,lwd=10*wm)
for(i in 1:5/10) g(i)
legend('center',bty='n',lwd=1
798::5,col=1:5, legend=c('10%','20%','30%','40%','50%')) pass=100 wm=0.2 curve( (x/(1-wm))/((pass-x)/wm),0,pass) uniroot(function(x,u0=1.2) x*wm/(1-wm)/(pass-x)-u0, c(0,pass)) (83/800)/(17/200) prop.test(c(83,17),c(800,200))
799:卵の名無しさん
18/10/14 07:14:30.48 +n7lCedi.net
pass=100
wm=0.2
curve( (x/(1-wm))/((pass-x)/wm),0,pass)
uniroot(function(x,u0=1.2) x*wm/(1-wm)/(pass-x)-u0, c(0,pass))
(83/800)/(17/200)
prop.test(c(83,17),c(800,200))
m_suc = function(rate) (pass* rate * (wm - 1))/(rate* (wm - 1) - wm)
m_suc(1.2)
curve(m_suc(x),1,5)
total=1000
wm=0.5
h <- function(mf){
M=ceiling(m_suc(mf))
F=pass-M
prop.test(c(M,F),c(total*(1-wm),total*wm))$p.value
}
r=seq(1,3,length=1000)
plot(r,sapply(r,h),type='l')
abline(h=0.05,lty=3)
uniroot(function(x)h(x)-0.05,c(1,2))$root
m_suc(1.43)
59/41
prop.test(c(59,41),c(500,500),correct=F)
prop.test(c(60,40),c(500,500),correct=F)
800:卵の名無しさん
18/10/17 18:14:18.05 Bhck6Le9.net
,イ ヽ
/ ヽ ヽ ヽ
/ / | {. ハ ハ ヽ
/ :! │ j 八 | | )
,′ i | ｛. i | {＼ ﾊ | i| ﾘ )､ / 国試浪人事務員は薄汚いジジイ
ｉ i | 八 ｉ Vx炒¨＾jﾉ}八/ ( ヽ / 敗者は敗者らしく
| ｉ 从 { { ＼八{ //// ｲ )ヾﾉ ( この世界のなるべく闇いところで
Viﾊ ＼|x炒＾j ( (ｲ / 謝りながら生きていくといいのよ
＼＼_V//く _ ) ) ヾ ゝ,
`ｰ小. _-ｧ ( ( ヾ
|i∧ く:::_ノ ) i 丶､ﾞゝ
|ｉ| {丶、 ´ ( | iヾﾞ
ﾉ八{ { {＞ . ／) ｉ ヾ
〃{｛ `ｰ=ミ ｀T爪八 i , ヾ 八
ヾ人＼ ｀ヽ j{__j ｼ｀ヾﾐﾐﾙ'ﾊヽ.ト＜7＼＿＿＿___
／⌒￢=-ミ#ヾ/ ＼ ＼/####/ ｀ヽ
/ ヽ##｀ヽ､ノ⌒ヽ V⌒ゝ､#;;| '.
{ ﾞi####ヽ 八 |'⌒´##ヾ |
| |#####＼/:i ＼'ー‐'i"######ﾞ) |
│ |########| :| ( ⌒ )#######| |
801:卵の名無しさん
18/10/21 21:20:01.39 Z8EDA9Rx.net
# 縦3マス、横4マスの12マスのうちランダムに選ばれた2マスにそれぞれ宝が眠っている。
# AEIBFJ…の順で縦に宝を探していく方法をとるP君と、
# ABCDEFGH…の順で横に宝を探していく方法をとるQ君が、同時に地点Aから探索を開始した。
# どっちの方が有利？
#
# ABCD
# EFGH
# IJKL
y=LETTERS[1:12]
z=matrix(y,ncol=4,byrow=T)
as.vector(z) # P
as.vector(t(z)) # Q(=y)
x=c(1,1,rep(0,10))
is.P1st <- function(){
Q=sample(x)
z=matrix(Q,ncol=4,byrow=T)
P=as.vector(z)
which.max(P) < which.max(Q)
}
k=1e3
re=replicate(k,mean(replicate(k,is.P1st())))
summary(re)
802:卵の名無しさん
18/10/21 21:45:36.09 Z8EDA9Rx.net
> x=c(1,1,rep(0,10))
> PQ <- function(){
+ Q=sample(x)
+ z=matrix(Q,ncol=4,byrow=T)
+ P=as.vector(z)
+ c( even=which.max(P) == which.max(Q),
+ p1st=which.max(P) < which.max(Q),
+ q1st=which.max(P) > which.max(Q))
+
+ }
> k=1e6
> re=replicate(k,PQ())
> mean(re['even',]) ; 13/(26+27+13)
[1] 0.197025
[1] 0.1969697
> mean(re['p1st',]) ; 26/(26+27+13)
[1] 0.393803
[1] 0.3939394
> mean(re['q1st',]) ; 27/(26+27+13)
[1] 0.409172
[1] 0.4090909
803:卵の名無しさん
18/10/22 02:17:14.19 YJqFFfdN.net
>>748
一般解
treasure <- function(m=3,n=4,k=2){
y=1:(m*n)
(z=matrix(y,ncol=n,byrow=T))
(P=as.vector(z))
(Q=as.vector(t(z)))
PQ <- function(x){
p=q=numeric(k)
for(i in 1:k){
p[i]=which(P==x[i])
q[i]=which(Q==x[i])
}
min(p)-min(q)
}
tre=combn(m*n,k)
re=apply(tre,2,PQ)
c(P1st=sum(re<0),Q1st=sum(re>0),even=sum(re==0))
}
treasure(3,4,2)
treasure(5,10,3)