###################################### # Statistical methods for data science: R # # http://didawiki.di.unipi.it/doku.php/mds/smd/ ###################################### # [CTRL-Enter] to execute a line ###### Script R12: Estimators and MSE ########## # compute pair of estimators (S, T) st.est = function(data) { s = sum(data==0)/length(data) t = exp(-mean(data)) return( c(s, t) ) } n = 100 mu = log(10) p0 = exp(-mu) arrivals = rpois(n, mu) arrivals data = replicate(1000, st.est(rpois(n, mu))) plot(density(data[2,]), col="blue", main='') lines(density(data[1,]), col="red") legend('topright', pch=16, c('S', 'T'), col=c('red', 'blue')) # textbook: Fig. 19.4 et = function(n, mu) exp(-n*mu*(1-exp(-1/n))) # see Exercise 19.9 plot(1:n, et(1:n, mu), ylim=c(0, .2), pch=16) abline(h=p0, col=2) # E[T]->p0 for n->infty # Var(S) = p0*(1-p0)/n where p0 = exp(-mu) var.s = function(n, mu) exp(-mu)*(1-exp(-mu))/n plot(1:n, var.s(1:n, mu), ylim=c(0, .1), pch=16) # Var(T) = E[T^2] - E[T]^2 var.t = function(n, mu) exp(-n*mu*(1-exp(-2/n))) - et(n, mu)^2 points(1:n, var.t(1:n, mu), ylim=c(0, 0.25), pch=16, col="red") # [T] textbook: Fig. 20.3 mse.s = function(n, mu) var.s(n, mu) + 0 mse.t = function(n, mu) var.t(n, mu) + (et(n, mu) - exp(-mu))^2 mus = seq(0, 5, 0.05) plot(mus, mse.s(n, mus), pch=16, cex=0.5, col="red", xlab=expression(mu), ylab='MSE') points(mus, mse.t(n, mus), pch=16, cex=0.5, col="blue") legend('topright', pch=16, c('S', 'T'), col=c('red', 'blue')) # The German Tank exercise t12.est = function(data) { n = length(data) t1 = 2*mean(data)-1 t2 = (n+1)/n*max(data)-1 #print(c(t1, t2)) return( c(t1, t2) ) } N = 1000 # tanks produced n = 100 # tanks destroyed data = replicate(2000, t12.est(sample(1:N, n))) # [T] textbook: Fig. 20.1 plot(density(data[2,]), col="red", main="Tanks' estimators", xlim=c(N*0.8, N*1.2)) lines(density(data[1,]), col="blue") legend('topleft', pch=16, c('T2', 'T1'), col=c('red', 'blue'))