Geometric distribution, sample size n, probability of success p.

To test whether the P_hat_mle is biased or not.

p_hat_mle = n/(n+sum(x_i)) = (total # of success)/(total # of trials)

p = seq(0,1,0.01)
p = p[2:100]
Ed1 = p*log(p)/(p-1)  # n=1, MLE
EdBayes = p*(3 - 4*p + p^2 + 2 *log(p))/(-1 + p)^3 #n=1 Bayes prior: beta(1,1)
Ed2 = (2*p*(1 - p + p* log(p)))/(-1 + p)^2  # n=2, MLE
#n=5, MLE
Ed5 = (5*p*(-3+16*p-36*p^2+48*p^3-25*p^4+12*p^4*log(p)))/(12*(-1+p)^5)
# n=10, MLE
Ed10 = (1/(252*(p - 1)^10))*(p*(-(7129*p^9) + 2520*p^9*log(p) + 22680*p^8 - 45360*p^7 + 70560*p^6 - 79380*p^5 + 63504*p^4 - 35280*p^3 + 12960*p^2 - 2835*p + 280))
#n=100
Ed100 = (1/ p)*(100*(1 -p)*( p^2/(99*(1 -p)^2) -p^3/(98*(1 -p)^3) +p^4/(97*(1 -p)^4) - p^5/(96*(1 -p)^5) +p^6/(95*(1 -p)^6) -p^7/(94*(1 -p)^7) +p^8/(93*(1 -p)^8) - p^9/(92*(1 -p)^9) +p^10/(91*(1 -p)^10) -p^11/(90*(1 -p)^11) + p^12/(89*(1 -p)^12) -p^13/(88*(1 -p)^13) +p^14/(87*(1 -p)^14) - p^15/(86*(1 -p)^15) +p^16/(85*(1 -p)^16) -p^17/(84*(1 -p)^17) + p^18/(83*(1 -p)^18) -p^19/(82*(1 -p)^19) +p^20/(81*(1 -p)^20) - p^21/(80*(1 -p)^21) +p^22/(79*(1 -p)^22) -p^23/(78*(1 -p)^23) + p^24/(77*(1 -p)^24) -p^25/(76*(1 -p)^25) +p^26/(75*(1 -p)^26) - p^27/(74*(1 -p)^27) +p^28/(73*(1 -p)^28) -p^29/(72*(1 -p)^29) + p^30/(71*(1 -p)^30) -p^31/(70*(1 -p)^31) +p^32/(69*(1 -p)^32) - p^33/(68*(1 -p)^33) +p^34/(67*(1 -p)^34) -p^35/(66*(1 -p)^35) + p^36/(65*(1 -p)^36) -p^37/(64*(1 -p)^37) +p^38/(63*(1 -p)^38) - p^39/(62*(1 -p)^39) +p^40/(61*(1 -p)^40) -p^41/(60*(1 -p)^41) + p^42/(59*(1 -p)^42) -p^43/(58*(1 -p)^43) +p^44/(57*(1 -p)^44) - p^45/(56*(1 -p)^45) +p^46/(55*(1 -p)^46) -p^47/(54*(1 -p)^47) + p^48/(53*(1 -p)^48) -p^49/(52*(1 -p)^49) +p^50/(51*(1 -p)^50) - p^51/(50*(1 -p)^51) +p^52/(49*(1 -p)^52) -p^53/(48*(1 -p)^53) + p^54/(47*(1 -p)^54) -p^55/(46*(1 -p)^55) +p^56/(45*(1 -p)^56) - p^57/(44*(1 -p)^57) +p^58/(43*(1 -p)^58) -p^59/(42*(1 -p)^59) + p^60/(41*(1 -p)^60) -p^61/(40*(1 -p)^61) +p^62/(39*(1 -p)^62) - p^63/(38*(1 -p)^63) +p^64/(37*(1 -p)^64) -p^65/(36*(1 -p)^65) + p^66/(35*(1 -p)^66) -p^67/(34*(1 -p)^67) +p^68/(33*(1 -p)^68) - p^69/(32*(1 -p)^69) +p^70/(31*(1 -p)^70) -p^71/(30*(1 -p)^71) + p^72/(29*(1 -p)^72) -p^73/(28*(1 -p)^73) +p^74/(27*(1 -p)^74) - p^75/(26*(1 -p)^75) +p^76/(25*(1 -p)^76) -p^77/(24*(1 -p)^77) + p^78/(23*(1 -p)^78) -p^79/(22*(1 -p)^79) +p^80/(21*(1 -p)^80) - p^81/(20*(1 -p)^81) +p^82/(19*(1 -p)^82) -p^83/(18*(1 -p)^83) + p^84/(17*(1 -p)^84) -p^85/(16*(1 -p)^85) +p^86/(15*(1 -p)^86) - p^87/(14*(1 -p)^87) +p^88/(13*(1 -p)^88) -p^89/(12*(1 -p)^89) + p^90/(11*(1 -p)^90) -p^91/(10*(1 -p)^91) +p^92/(9*(1 -p)^92) - p^93/(8*(1 -p)^93) +p^94/(7*(1 -p)^94) -p^95/(6*(1 -p)^95) + p^96/(5*(1 -p)^96) -p^97/(4*(1 -p)^97) +p^98/(3*(1 -p)^98) - p^99/(2*(1 -p)^99) +p^100/(1 -p)^100 + ( p^101*log(p))/(1 -p)^101))


par(mar=c(5,3,1,1)+0.1)
plot(p,Ed1,type='l',col='red', asp=1) #n=1
lines(p,p,type='l', col='pink')
lines(p,EdBayes,type='l', col='blue') #bayes, beta(1,1)
lines(p,Ed2,type='l', col='green') #n=2
lines(p,Ed5,type='l', col='purple') #n=5
lines(p,Ed10,type='l', col='orange') #n=10
lines(p,Ed100,type='l', col='black') #n=100

# This figure also shows a hint of information I(p)
# I suspect the mess when p->1 is round off error of R

plot(p,Ed1-p, type='l',ylim=c(-0.4,0.4),col='red')
abline(0,0,col='pink')
lines(p,EdBayes-p,type='l', col='blue')
lines(p,Ed2-p,type='l', col='green')
lines(p,Ed5-p,type='l', col='purple')
lines(p,Ed10-p,type='l', col='orange')
lines(p,Ed100-p,type='l', col='black')

 (Ed100-p)[1:50] # Ed100 is very close to p! ---- asymptotic unbiased
##  [1] 0.0000999796 0.0001978990 0.0002937596 0.0003875625 0.0004793091
##  [6] 0.0005690007 0.0006566384 0.0007422236 0.0008257575 0.0009072413
## [11] 0.0009866765 0.0010640641 0.0011394054 0.0012127018 0.0012839545
## [16] 0.0013531646 0.0014203336 0.0014854625 0.0015485528 0.0016096055
## [21] 0.0016686220 0.0017256035 0.0017805512 0.0018334664 0.0018843504
## [26] 0.0019332043 0.0019800294 0.0020248269 0.0020675981 0.0021083441
## [31] 0.0021470663 0.0021837658 0.0022184440 0.0022511019 0.0022817408
## [36] 0.0023103620 0.0023369666 0.0023615559 0.0023841312 0.0024046935
## [41] 0.0024232442 0.0024397844 0.0024543154 0.0024668384 0.0024773545
## [46] 0.0024858650 0.0024923711 0.0024968740 0.0024993749 0.0024998750