# example 3.22 in Davison(2008)
r = 2; m = 10; mu = 0; sig2 = 1;
return(r * y - m * log(1 + exp(y)) - (y - mu)^2 / (2 * sig2))
return(r * y - m * log.(1 .+ exp.(y)) .- (y .- mu) .^2 ./ (2 * sig2))
return(r - m * exp(y) / (1 + exp(y)) - (y - mu) / sig2)
return(r .- m * exp.(y) ./ (1 .+ exp.(y)) .- (y .- mu)./sig2)
function zfix(yfixed::Array)
zfixed = yf0 .+ (h(yf0) .- h(yf1) .+ (yf1 .- yf0) .* dh(yf1)) ./ (dh(yf1) .- dh(yf0))
# evaluate log-density (not necessary)
function hplus(y::Float64, yfixed::Array)
if i == 1 && y < zfixed[i]
return(h(yfixed[i]) + (y - yfixed[i]) * dh(yfixed[i]))
elseif i < n && y >= zfixed[i] && y < zfixed[i+1]
return(h(yfixed[i+1]) + (y - yfixed[i+1]) * dh(yfixed[i+1]))
elseif i == n && y >= zfixed[n]
return(h(yfixed[i]) + (y - yfixed[i]) * dh(yfixed[i]))
function gplus_cdf(yfixed::Array, zfixed::Array)
## integral from -infty to zi
# s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - )
s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - exp((yl-yfixed[i]) * dh(yfixed[i])))
s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((yu-yfixed[i]) * dh(yfixed[i])) - exp((zfixed[n]-yfixed[i]) * dh(yfixed[i])))
s[i] = exp(h(yfixed[i])) / dh(yfixed[i]) * (exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])) - exp((zfixed[i]-yfixed[i]) * dh(yfixed[i])))
return cumsum(pr), sum(s)
# sample from gplus density
function gplus_sample(yfixed)
gp = gplus_cdf(yfixed, zfixed)
ey = u * dh(yfixed[i]) * norm_const / exp(h(yfixed[i])) + exp((yl-yfixed[i])*dh(yfixed[i]))
return(yfixed[i] + log(ey)/dh(yfixed[i]))
elseif i == n && u >= zpr[i]
ey = (u - zpr[i]) * dh(yfixed[i+1]) * norm_const / exp(h(yfixed[i+1])) + exp((zfixed[i]-yfixed[i+1])*dh(yfixed[i+1]))
return(yfixed[i+1] + log(ey)/dh(yfixed[i+1]))
elseif i < n && u >= zpr[i] && u < zpr[i+1]
ey = (u - zpr[i]) * dh(yfixed[i+1]) * norm_const / exp(h(yfixed[i+1])) + exp((zfixed[i]-yfixed[i+1])*dh(yfixed[i+1]))
return(yfixed[i+1] + log(ey)/dh(yfixed[i+1]))