return(b * sum(x .* y) - sum(log.(1 .+ exp.(b*x))) - sum(1 .- 1 ./ (1 .+ exp.(b*x))) - b^2/(2*tau^2))
return(sum(x .* y) - sum(x .* (1 .- 1 ./ (1 .+ exp.(b*x)) )) + sum(x .* exp.(b*x) ./ (1 .+ exp.(b*x)).^2 ) - b/tau^2)
function phi(x, yfixed::Array)
if exp(h(x)) <= exp(hplus(x, yfixed))
return(exp(hplus(x, yfixed)))
x = arms(10000, [-1.5, -0.7, 0.7, 1.5])