expt.partly_comp1.r

## harmonic number ##
harmonic <- function(n){
  sum(sapply(seq(1,n),function(n){ return(1/n) }))
}
## expectation of partly complete case ##
expt.partly_comp1 <- function(m,n){
  if (n < m) {  
    m * ( harmonic(m) - harmonic(m-n) )
  }else{
    m * harmonic(m)
  }
}
## AKB48 ポスター部分コンプリート例 ##
> expt.partly_comp1(44,44)
# [1] 192.3999
> expt.partly_comp1(44,10)
# [1] 11.1987
> expt.partly_comp1(44,5)
# [1] 5.244046
## case when m = 20, n in [1,20] ##
seq.expt.partly_comp1 <- function(m){
    s <- seq(1, m)
    expt.partly_comp1 <- function(n){
        if (n < m) {  
            m * ( harmonic(m) - harmonic(m-n) )
        }else{
            m * harmonic(m)
        }
    }
    sapply(s, expt.partly_comp1)
}
seq.expt.partly_comp1(20)
# [1]    1.000000  2.052632  3.163743  4.340213  5.590213
# [5]    6.923547  8.352118  9.890580 11.557246 13.375428
# [10]  15.375428 17.597650 20.097650 22.954793 26.288126
# [15]  30.288126 35.288126 41.954793 51.954793 71.954793
## plot ##
plot(seq.expt.partly_comp1(20),type='l',xlab='n',ylab='expt.times')
plot(seq.expt.partly_comp1(100),type='l',xlab='n',ylab='expt.times')