Density 2d from curves

Question

I have a multi-parameter function on which I output parameters using MCMC. This means that I have many sample parameters, and I can build functions:

# Simulate some parameters. Really, I get these from MCMC sampling.
first = rnorm(1000)  # a
second = rnorm(1000)  # b

# The function (geometric)
geometric = function(x, a, b) b*(1 - a^(x + 1)/a)  

# Plot curves. Perhaps not the most efficient way, but it works.
curve(geometric(x, first[1], second[1]), ylim=c(-3, 3))  # first curve
for(i in 2:length(first)) {
  curve(geometric(x, first[i], second[i]), add=T, col='#00000030')  # add others
}

How can I do this on a density graph, and not on a graph of individual curves? For example, it is difficult to understand how denser it is y = 0 than around other values.

It would be nice:

• Ability to draw observable values ​​from above (points and lines).
• Drawing a contour line in density, for example. 95% highest attenuation or 2.5 or 97.5 quantiles.