How to describe/explain the shape of a distribution which has two peaks? How can I describe the shape of the following distribution, which has two peaks? What are the important things in describing the distribution shape?

It is the output of the following R code:
plot(density(Boston$tax))
(The data set is the Boston data from the MASS package.)
I would also like to know what the rug (these red points) means in the following plot:

I created it with the following code using same data as above:
plot(density(Boston$tax))
rug(Boston$tax, col=2, lwd=3.5)

 A: To describe such a two-peak shape to another person, you'd call it 'bimodal' (which just means 'two modes' - generally taken to be two local modes, even though only one of them might be 'the mode' of the distribution).
You could then seek to describe the locations and spreads and relative proportions or heights of the peaks (this might be done visually, or more formally, for example with Gaussian mixture model).
e.g. as a first, simple description, I might say "the distribution appears to be bimodal with the main peak at around 290 and a lower peak around 670"
- and then if necessary give additional detail on the relative heights at or widths-of/ areas-under the curve around the peaks, if any of those details matter to your audience (e.g. something along the lines of "the spread of the peak around the main mode is wider than that around the smaller mode").
If we're thinking in terms of something close to a Gaussian mixture model, there's a slight suggestion of a third "bump" coming in near 420, but its close enough to the bigger mode that it doesn't make a separate peak.

Those red marks you got using rug are the actual data values; for each observation, one red mark is placed in the margin (akin to the marks you see with stripchart(Boston$tax,pch="|")). You'd normally mark them with thin lines rather than wide lines as you have there. Because the values are placed at the edge of a plot, the marks look a little like the fringe tassels on the edges of a rug. That has nothing to do with the kernel density estimate itself (other than showing the data from which the KDE was computed), it's just adding a different kind of information to the plot; you could use rug to add information to various other displays of the data. A rug-plot is just a marginal (to some main plot) one-dimensional plot of the data values.
A rug plot is not as informative when the data has a lot of repeated values (as the tax variable does - for example, there are 132 values at 666 which are all drawn one on top of the other); you can normally improve that with a little bit of 'jitter', but there are so many repeated values that even rug(jitter(x,amount=20)) doesn't distinguish the values. For that situation might be better to use transparency and a smaller amount of jitter, or some other indication, such as a dot plot:
plot(density(Boston$tax),col=3)  #$
stripchart(Boston$tax,add=TRUE,pch=16,cex=.03,at=-0.00013,method="stack",col=8)


