expt.partly_comp2.r

as.prob <- function(weight){
  return(weight/sum(weight))
}
expt.partly_comp2 <- function(vec.prob, n){
  m <- length(vec.prob)
  g <- function(m, n, i, vec.prob){
    (-1)**(n-1-i)*choose(m-i-1, m-n)*1/(1-apply(combn(vec.prob, i), 2, sum))
  }
  tmp <- sapply(seq(0,n-1),function(i){
    g( m, n, i, vec.prob)
  })
  res <- sum(sapply(tmp,sum))
  return(res)
} 

## 等確率
weight <- rep(1, 20)
expt.partly_comp2(as.prob(weight),20) # 完全コンプの場合
# [1] 71.95479 ※ 全節と同じ結果

## レアが10倍出にくい（レア：10個，通常：10個）
weight <- c(rep(1, 10),rep(10, 10))
expt.partly_comp2(as.prob(weight),20)  # 完全コンプの場合
# [1] 322.1871

## レアが100倍出にくい（レア：10個，通常：10個）
weight <- c(rep(1, 10),rep(100, 10))
expt.partly_comp2(as.prob(weight),20)
# [1] 2958.258

seq.expt.partly_comp2 <- function(vec.prob){
    m <- length(vec.prob)
    expt.partly_comp2 <- function(n){
       m <- length(vec.prob)
       g <- function(m, i, vec.prob){
          (-1)**(n-1-i)*choose(m-i-1, m-n)*1/(1-apply(combn(vec.prob, i), 2, sum))
       }
       tmp <- sapply(seq(0,n-1),function(i){
         g( m, i, vec.prob)
       })
       res <- sum(sapply(tmp,sum))
       return(res)
    }
    sapply(seq(1, m), expt.partly_comp2)
}

partial = c(8,9,10,11,12)

# レアが10倍出にくい（レア：10個，通常：10個）
weight <- c(rep(1, 10),rep(10, 10))
s1 <- seq.expt.partly_comp2(as.prob(weight))
s1
# [1]  1.000   2.091   3.293   4.627   6.126
# [5]  7.834   9.816  12.168  15.053  18.761  
#[10] 23.906  31.266  41.474  54.932  72.008
#[15] 93.397 120.641 157.217 212.192 322.187
plot(s1, type='b',xlab='b',ylab='expt.times')
plot(partial, s1[partial], type='b',xlab='n',ylab='expt.times')

# レアが50倍出にくい（レア：10個，通常：10個）
weight <- c(rep(1, 10),rep(50, 10))
s2 <- seq.expt.partly_comp2(as.prob(weight))
s2
# [1]    1.000    2.106    3.345    4.750    6.375
# [5]    8.298   10.653   13.690   17.961   25.215 
#[10]   57.885  109.104  171.658  244.307  329.277  
#[15]  431.274  558.773  728.773  983.773 1493.773
plot(s2, type='b',xlab='b',ylab='expt.times')
plot(partial, s2[partial], type='b',xlab='n',ylab='expt.times')

# レアが100倍出にくい（レア：10個，通常：10個）
weight <- c(rep(1, 10),rep(100, 10))
s3 <- seq.expt.partly_comp2(as.prob(weight))
s3
# [1]    1.000    2.108    3.352    4.769    6.414
# [5]    8.374   10.798   13.973   18.575   26.989 
#[10]  105.131  213.687  339.513  483.760  652.091
#[15]  854.091 1106.591 1443.257 1948.257 2958.257
plot(s3, type='b',xlab='b',ylab='expt.times')
plot(partial, s3[partial], type='b',xlab='n',ylab='expt.times')

partial = c(3,4,5,6,7)

# レアが10倍出にくい（レア：15個，通常：5個）
weight <- c(rep(1, 15),rep(10, 5))
s4 <- seq.expt.partly_comp2(as.prob(weight))
s4
# [1]   1.000   2.143   3.473   5.052   6.977   9.398  12.454
# [8]  16.235  20.781  26.108  32.243  39.251  47.258  56.483
#[15]  67.287  80.275  96.520 118.185 150.684 215.684
plot(s4, type='b',xlab='b',ylab='expt.times')
plot(partial, s4[partial], type='b',xlab='n',ylab='expt.times')

# レアが50倍出にくい（レア：15個，通常：5個）
weight <- c(rep(1, 15),rep(50, 5))
s5 <- seq.expt.partly_comp2(as.prob(weight))
s5
# [1]   1.000   2.219   3.779   5.941   9.462  20.708  37.410
# [8]  57.170  79.103 103.161 129.655 159.098 192.223 230.080
#[15] 274.247 327.247 393.497 481.830 614.330 879.330
plot(s5, type='b',xlab='n',ylab='expt.times')
plot(partial, s5[partial], type='b',xlab='n',ylab='expt.times')

# レアが100倍出にくい（レア：15個，通常：5個）
weight <- c(rep(1, 15),rep(100, 5))
s6 <- seq.expt.partly_comp2(as.prob(weight))
s6
# [1]    1.000    2.234    3.843    6.155   10.273   36.246
# [7]   71.406  110.771  153.655  200.469  251.969  309.191
#[13]  373.566  447.137  532.971  635.971  764.721  936.387
#[19] 1193.887 1708.887
plot(s6, type='b',xlab='n',ylab='expt.times')
plot(partial, s6[partial], type='b',xlab='n',ylab='expt.times')