mean n binom

###########
### comparison of simulated cumulative sum of
### geometric variables with mean = 0
### (ie. cumsum(geomVariable(p)-1/p))
### with a mean of negative binomial distributions
###vars definition
D        <- 1e06                      # D=1/p
y        <- c(2, 5, 10, 15, 20, 25)   # n to test
x.length <- 200                       # number of x's to test
x<-seq(-D*10,D*10,length=x.length)    # sequence of x.length variables from min and max x
                                      # x being the the cumulative sum assessed

### dnb.mean.func calculates the mean negative binomial density
### for 1:n sucesses after making mean = 0
###
dnb.mean.func<-function(x,y,D){
    z<-matrix(0,y,1)
    for(i in 1:y){
    z[i]<-dnbinom(floor(x+i*(1-1/D)/(1/D)),i,1/D)
    }
    return(mean(z))
}

plot.3.1.mat<-matrix(0, x.length, 3);plot.3.1.mat[,1]<-x
plot.3.2.mat<-matrix(0, x.length ,3);plot.3.2.mat[,1]<-x
plot.3.3.mat<-matrix(0, x.length ,3);plot.3.3.mat[,1]<-x
plot.3.4.mat<-matrix(0, x.length, 3);plot.3.4.mat[,1]<-x
plot.3.5.mat<-matrix(0, x.length ,3);plot.3.5.mat[,1]<-x
plot.3.6.mat<-matrix(0, x.length ,3);plot.3.6.mat[,1]<-x

###assessing the mean negative binomial density at x

for(i in 1:200){
    plot.3.1.mat[i,2]<-dnb.mean.func(plot.3.1.mat[i,1],y[1],D) #density values for first n to test
    plot.3.2.mat[i,2]<-dnb.mean.func(plot.3.2.mat[i,1],y[2],D) #density values for second n to test
    plot.3.3.mat[i,2]<-dnb.mean.func(plot.3.3.mat[i,1],y[3],D) #density values for third n to test
    plot.3.4.mat[i,2]<-dnb.mean.func(plot.3.4.mat[i,1],y[4],D) #density values for fourth n to test
    plot.3.5.mat[i,2]<-dnb.mean.func(plot.3.5.mat[i,1],y[5],D) #density values for fifth n to test
    plot.3.6.mat[i,2]<-dnb.mean.func(plot.3.6.mat[i,1],y[6],D) #density values for sixth n to test
    print(i)
}

###simulation of cumulative sum - D
###creates many simulations of a cumulative sum for each n
###can be normalised to btc instead of number of fails (shares submitted)
###
iter   <- 5e05      #number of iterations of the cumulative sum,
                    #more = less variability in result, takes longer
rgIter <- vector('integer', (iter*3))
rg.mat <- c(rep(y[1], iter),rep(y[2], iter),rep(y[3], iter),
            rep(y[4], iter),rep(y[5], iter),rep(y[6], iter))

rg<-function(n){cumsum(rgeom(n,1/D)-D)}

###applies rg to rg.mat, taking n from rg.mat[,1]
rgIter <- mapply(rg,rg.mat)

###creates vectors from rgIter list
rgIter1 <- unlist(rgIter[1:iter])
rgIter2 <- unlist(rgIter[(iter+1):(iter*2)])
rgIter3 <- unlist(rgIter[(iter*2+1):(iter*3)])
rgIter4 <- unlist(rgIter[(iter*3+1):(iter*4)])
rgIter5 <- unlist(rgIter[(iter*4+1):(iter*5)])
rgIter6 <- unlist(rgIter[(iter*5+1):(iter*6)])
rgIter<-0;rg.mat<-0

###collecting simulation density values for quicker plotting
rgDens1 <- density(rgIter1)
rgDens2 <- density(rgIter2)
rgDens3 <- density(rgIter3)
rgDens4 <- density(rgIter4)
rgDens5 <- density(rgIter5)
rgDens6 <- density(rgIter6)

###simple plot of simulation densities and 'mean' negative binomial densities
png('mean-dnb-plot.png', height=640*3,width=640*16/9,res=130)
par(mfrow=c(3, 2))

plot(plot.3.1.mat[,1], plot.3.1.mat[, 2],type='l',col='red',lwd=1.5,xlim=range(x),xlab='',ylab='')
lines(rgDens1$x,rgDens1$y,type='l',lwd=1)
legend("topright",legend=c('Simulation\ndensity', 'Mean\nnegative\nbinomial\nprobability\ndensities'),col=c('black','red'),lwd=c(1,1),bty='n')
title(main=paste(y[1],"sucesses"),xlab="Cumulative sum of total round shares minus D", ylab="Density")

plot(plot.3.2.mat[,1],plot.3.2.mat[,2],type='l',col='red',lwd=1.5,xlim=range(x),xlab='',ylab='')
lines(rgDens2$x,rgDens2$y,type='l',lwd=1)
legend("topright",legend=c('Simulation\ndensity', 'Mean\nnegative\nbinomial\nprobability\ndensities'),col=c('black','red'),lwd=c(1,1),bty='n')
title(main=paste(y[2],"sucesses"),xlab="Cumulative sum of total round shares minus D", ylab="Density")

plot(plot.3.3.mat[,1],plot.3.3.mat[,2],type='l',col='red',lwd=1.5,xlim=range(x),xlab='',ylab='')
lines(rgDens3$x,rgDens3$y,type='l',lwd=1)
legend("topright",legend=c('Simulation\ndensity', 'Mean\nnegative\nbinomial\nprobability\ndensities'),col=c('black','red'),lwd=c(1,1),bty='n')
title(main=paste(y[3],"sucesses"),xlab="Cumulative sum of total round shares minus D", ylab="Density")

plot(plot.3.4.mat[,1],plot.3.4.mat[,2],type='l',col='red',lwd=1.5,xlim=range(x),xlab='',ylab='')
lines(rgDens4$x,rgDens4$y,type='l',lwd=1)
legend("topright",legend=c('Simulation\ndensity', 'Mean\nnegative\nbinomial\nprobability\ndensities'),col=c('black','red'),lwd=c(1,1),bty='n')
title(main=paste(y[4],"sucesses"),xlab="Cumulative sum of total round shares minus D", ylab="Density")

plot(plot.3.5.mat[,1],plot.3.5.mat[,2],type='l',col='red',lwd=1.5,xlim=range(x),xlab='',ylab='')
lines(rgDens5$x,rgDens5$y,type='l',lwd=1)
legend("topright",legend=c('Simulation\ndensity', 'Mean\nnegative\nbinomial\nprobability\ndensities'),col=c('black','red'),lwd=c(1,1),bty='n')
title(main=paste(y[5],"sucesses"),xlab="Cumulative sum of total round shares minus D", ylab="Density")

plot(plot.3.6.mat[,1],plot.3.6.mat[,2],type='l',col='red',lwd=1.5,xlim=range(x),xlab='',ylab='')
lines(rgDens6$x,rgDens6$y,type='l',lwd=1)
legend("topright",legend=c('Simulation\ndensity', 'Mean\nnegative\nbinomial\nprobability\ndensities'),col=c('black','red'),lwd=c(1,1),bty='n')
title(main=paste(y[6],"sucesses"),xlab="Cumulative sum of total round shares minus D", ylab="Density")

dev.off()