#	Example 2:

m 	<- 1000
theta	<- 0.5
alpha	<- 0.05
p	<- 0.05
#p	<- 0.05/1000

theta.vec	<- rep(theta,m)
t.vec		<- sort(rnorm(m,theta,1))

sel.vec	<- qnorm(1 - alpha/2) < abs(t.vec)
cov.vec	<- abs(t.vec - theta.vec) < qnorm(1-p/2)

plt.green	<- sel.vec & (cov.vec)
plt.red	<- sel.vec & (!cov.vec)

plot(t.vec,ylim=c(theta-4,theta+4),main = paste(" R = ",sum(sel.vec),", V = ",sum(plt.red)) )

segments((1:m)[plt.red],(t.vec - qnorm(1-p/2))[plt.red],(1:m)[plt.red],(t.vec + qnorm(1-p/2))[plt.red],col=2)
segments((1:m)[plt.green],(t.vec - qnorm(1-p/2))[plt.green],(1:m)[plt.green],(t.vec + qnorm(1-p/2))[plt.green],col=3)
abline(h = theta, col=4,lwd = 2)




#	Example 3:


m 		<- 100
#m 		<- 10000
pi.1		<- 0.8
pi.2		<- 0.2
theta.1	<- -0.2
theta.2	<- 4

theta.vec	<- c(rep(theta.1,pi.1*m),rep(theta.2,pi.2*m))
t.vec		<- rnorm(m,theta.vec)

theta.srt	<- theta.vec[order(t.vec)]
t.srt		<- t.vec[order(t.vec)]
p.srt		<- 2*(1 - pnorm(abs(t.srt)))

#sel.vec	<- p.srt <= 0.05			# Marginal 0.05 selection
sel.vec	<- p.srt <= 0.05/(2*m)		# Bonferroni 0.05 selection
#sel.vec	<- p.adjust(p.srt,"BH") <= 0.05	# BH  0.05 selection

R		<- sum(sel.vec)
FCR.q		<- qnorm(1-(R*0.05)/(2*m))
cov.vec	<- abs(t.srt - theta.srt) < FCR.q
plt.green	<- sel.vec & (cov.vec)
plt.red	<- sel.vec & (!cov.vec)
V		<- sum(plt.red)
plot(t.srt,ylim=c(theta.1-4,theta.2+4),main = paste(" FCP = ",V,"/",R," = ",round(V/R,3)))
points(theta.srt,col=4)

segments((1:m)[plt.red],(t.srt -FCR.q)[plt.red],(1:m)[plt.red],(t.srt + FCR.q)[plt.red],col=2)
segments((1:m)[plt.green],(t.srt - FCR.q)[plt.green],(1:m)[plt.green],(t.srt + FCR.q)[plt.green],col=3)
abline(h = 0, col=4,lty=2)





#	Answer to PM's question

m 		<- 100
theta.vec	<- rexp(m,0.5)*sample(c(-1,1),m,replace = T)
t.vec		<- sort(rnorm(m,theta.vec))

#	STEP 1

sel.vec	<- rep(T,m)
FCR.q		<- qnorm(1-(sum(sel.vec)*0.05)/(2*m))
plt.green	<- sel.vec & FCR.q <= abs(t.vec)
plt.red	<- sel.vec &  FCR.q > abs(t.vec)

plot(t.vec,ylim=c(min(theta.vec)-4,max(theta.vec)+4),main = paste(" R = ",sum(sel.vec)))
segments((1:m)[plt.red],(t.vec -FCR.q)[plt.red],(1:m)[plt.red],(t.vec + FCR.q)[plt.red],col=2)
segments((1:m)[plt.green],(t.vec - FCR.q)[plt.green],(1:m)[plt.green],(t.vec + FCR.q)[plt.green],col=3)
abline(h = 0, col=1,lty=4)

plot(t.vec,ylim=c(min(theta.vec)-4,max(theta.vec)+4),main = paste(" R = ",sum(plt.green)))
segments((1:m)[plt.green],(t.vec - FCR.q)[plt.green],(1:m)[plt.green],(t.vec + FCR.q)[plt.green],col=3)
abline(h = 0, col=1,lty=4)


#	STEP 2,3, 4, .....

sel.vec	<- plt.green
FCR.q		<- qnorm(1-(sum(sel.vec)*0.05)/(2*m))
plt.green	<- sel.vec & FCR.q <= abs(t.vec)
plt.red	<- sel.vec &  FCR.q > abs(t.vec)

plot(t.vec,ylim=c(min(theta.vec)-4,max(theta.vec)+4),main = paste(" R = ",sum(sel.vec)))
segments((1:m)[plt.red],(t.vec -FCR.q)[plt.red],(1:m)[plt.red],(t.vec + FCR.q)[plt.red],col=2)
segments((1:m)[plt.green],(t.vec - FCR.q)[plt.green],(1:m)[plt.green],(t.vec + FCR.q)[plt.green],col=3)
abline(h = 0, col=1,lty=4)

#	(STOP WHEN sum(plt.red) = 0)






#	Example 4

q		<- 0.05
m 		<- 10^5

#theta.vec	<- c(rexp(0.9*m,1),rexp(0.1*m,.1))*sample(c(-1,1),m,replace = T)
#theta.vec	<- c(rep(0,0.9*m),rexp(0.1*m,.1))*sample(c(-1,1),m,replace = T)
theta.vec	<- c(rexp(0.9*m,3),rexp(0.1*m,1))*sample(c(-1,1),m,replace = T)
t.vec		<- rnorm(m,theta.vec)
p.vec		<- 2*(1 - pnorm(abs(t.vec)))
sel.vec	<- p.adjust(p.vec,"BH") <= q	

#	4.1 two sided CI's for all Theta's

R		<- m
plt.green	<- abs(t.vec - theta.vec) <= qnorm(1 - q/2)
plt.red	<- abs(t.vec - theta.vec) > qnorm(1 - q/2)
V		<- sum(plt.red)
plot(t.vec[plt.green],theta.vec[plt.green],col=3,ylab = "Theta", xlab = "T",
main = paste("CI's for all ",m," Theta's  -- FCP = ",V,"/",R," = ",round(V/R,4)))
points(t.vec[plt.red],theta.vec[plt.red],col=2)
abline(h =0,col=2,lty = 2)
abline(qnorm(1-q/2),1,col=4)
abline(-qnorm(1-q/2),1,col=4)


#	4.2 two sided CI's for BH selected Theta's

R		<- sum(sel.vec)
plt.green	<- sel.vec & (abs(t.vec - theta.vec) <= qnorm(1 - (R*q)/(2*m)))
plt.red	<- sel.vec & (abs(t.vec - theta.vec) > qnorm(1 - (R*q)/(2*m)))
V		<- sum(plt.red)
plot(t.vec[plt.green],theta.vec[plt.green],col=3,ylab = "Theta", xlab = "T",
main = paste("CI's for",R," BH selected Theta's  -- FCP = ",V,"/",R," = ",round(V/R,3)))
points(t.vec[plt.red],theta.vec[plt.red],col=2)
abline(h =0,col=2,lty = 2)
abline(qnorm(1 - (R*q)/(2*m)),1,col=4)
abline(-qnorm(1 - (R*q)/(2*m)),1,col=4)
abline(v = qnorm(1 - (R*q)/(2*m)),col=4,lty=2)
abline(v = -qnorm(1 - (R*q)/(2*m)),col=4,lty=2)



#	4.3 sign determining CI's for BH selected Theta's

plt.green	<- sel.vec & (sign(t.vec) == sign(theta.vec))
plt.red	<- sel.vec & (sign(t.vec) != sign(theta.vec))
V		<- sum(plt.red)
plot(t.vec[plt.green],theta.vec[plt.green],col=3,ylab = "Theta", xlab = "T",
main = paste("CI's for",R," BH selected Theta's  -- FCP = ",V,"/",R," = ",round(V/R,3)))
points(t.vec[plt.red],theta.vec[plt.red],col=2)
abline(h =0,col=2,lty = 2)
segments(qnorm(1 - (R*q)/(2*m)),0,40,0,col=4,lwd = 2)
segments(-qnorm(1 - (R*q)/(2*m)),0,-40,0,col=4,lwd = 2)
abline(v = qnorm(1 - (R*q)/(2*m)),col=4,lty=2)
abline(v = -qnorm(1 - (R*q)/(2*m)),col=4,lty=2)