# functions used to generate from conditional prior of nu given psi (need to specify z0 see notes)

# G0 cdf
G0=function(z0,z){
if (z < z0){
return(cat("argument to G0 wrong","\n"))}
G0=z0*(pnorm(z)-pnorm(z0))+(dnorm(z)-dnorm(z0))
G0=G0/(z0*(1-pnorm(z0))-dnorm(z0))
return(G0)
}

# G1 cdf
G1=function(z0,z){
if (z > z0){
return(cat("argument to G1 wrong","\n"))}
G1=z0*pnorm(z)+dnorm(z)
G1=G1/(z0*pnorm(z0)+dnorm(z0))
return(G1)
}


#---------------------------------------------------------------
# generate from density proportional to |z-z0|phi(z) (here w is a value from the uniform(0,1) distribution) 
# so nu=nu_{0)(psi) + tau_{20}(psi)*linnorm(z0,w) is the value from the conditional prior of nu given psi

linnorm=function(z0,w){

norm1=pnorm(z0,0,1)+dnorm(z0)/z0
norm0=-1+pnorm(z0,0,1)+dnorm(z0)/z0
norm=norm1+norm0
pz0=norm1/norm
ind=rbinom(1,1,pz0)

if (ind==0){
# lower bound zlow for bisection
zlow=z0
# find upper bound zup for bisection
zup=0
test=G0(z0,zup)
while (test < w ) {
zup=zup+abs(z0)
test=G0(z0,zup)
}
# now solve G0(z0,z)=w for z
z=(zlow+zup)/2
test=G0(z0,z)
while (abs(test-w)> 0.0001) {
if (test <= w){zlow=z
z=(zlow+zup)/2
test=G0(z0,z)}
if (test >= w){zup=z
z=(zlow+zup)/2
test=G0(z0,z)}
}
}

if (ind==1){
# upper bound zup for bisection
zup=z0
# find lower bound zlow for bisection
zlow=0
test=G1(z0,zlow)
while (test > w) {
zlow=zlow-abs(z0)
test=G1(z0,zup)
}
# now solve G1(z0,z)=w for z
z=(zlow+zup)/2
test=G1(z0,z)
while (abs(test-w)> 0.0001) {
if (test <= w){zlow=z
z=(zlow+zup)/2
test=G1(z0,z)}
if (test >= w){zup=z
z=(zlow+zup)/2
test=G1(z0,z)}
}
}
return(z)
}