fs = function(x,epsilon,delta) dnorm(sinh(delta*asinh(x)-epsilon))*delta*cosh(delta*asinh(x)-epsilon)/sqrt(1+x^2)
vec = seq(-5,5,0.001)
plot(vec,fs(vec,0,1.2),type="l",ylim=c(0,0.5),xlab="", ylab="Densité de probabilité", col="red")
points(vec,fs(vec,0,0.9),type="l")