n=10

 p<-seq(0,1,by=.01)
 msep = n*p*(1-p)/(n^2)
 msepstar = n*p*(1-p)/((n+2)^2) + ( (n*p+1)/(n+2) -p   )^2
 plot(p, msep-msepstar,type='l')


 n=100
 msep = n*p*(1-p)/(n^2)
 msepstar = n*p*(1-p)/((n+2)^2) + ( (n*p+1)/(n+2) -p   )^2
 plot(p, msep-msepstar,type='l')
 abline(h=0)