m  <- 1000
pi <-  rbeta(m,.5,.5)
n  <-  49
sqrtn <- sqrt(n)
pcut <- sqrt(2*sqrtn+1)/(2*(sqrtn+1))

x<- 1:m
for (i in 1:m) x[i]  <-  rbinom(1,n,pi[i])


phat <- x/n
ptilde <- (sqrtn*phat + .5)/(sqrtn + 1)

pest <- phat

padapt1 <- ifelse(abs(pest-.5)< pcut, ptilde, phat)
pest<-padapt1 
padapt2 <- ifelse(abs(pest-.5)< pcut, ptilde, phat)

gamma <- (.5-pest)^2/(pest*(1-pest)/n + (.5 -pest)^2)
pstar1 <- gamma*phat + .5*(1-gamma)
pest <- pstar1
gamma <- (.5-pest)^2/(pest*(1-pest)/n + (.5 -pest)^2)
pstar2 <- gamma*phat + .5*(1-gamma)


msephat <- mean((phat-pi)**2)
mseptilde <- mean((ptilde-pi)**2)
msepadapt1 <- mean((padapt1-pi)**2)
msepadapt2 <- mean((padapt2-pi)**2)
msepstar1 <- mean((pstar1-pi)**2)
msepstar2 <- mean((pstar2-pi)**2)
maephat <- mean(abs(phat-pi))
maeptilde <- mean(abs(ptilde-pi))
maepadapt1 <- mean(abs(padapt1-pi))
maepadapt2 <- mean(abs(padapt2-pi))
maepstar1 <- mean(abs(pstar1-pi))
maepstar2 <- mean(abs(pstar2-pi))

msephat
mseptilde
msepadapt1 
msepadapt2 
msepstar1
msepstar2

maephat
maeptilde
maepadapt1
maepadapt2 
maepstar1
maepstar2

var(x)