Can somebody offer an example of a unimodal distribution which has a skewness of zero but which is not symmetrical? In May 2010 Wikipedia user Mcorazao added a sentence to the skewness article that "A zero value indicates that the values are relatively evenly distributed on both sides of the mean, typically but not necessarily implying a symmetric distribution." However, the wiki page has no actual examples of distributions which break this rule. Googling "example asymmetrical distributions with zero skewness" also gives no real examples, at least in the first 20 results. 
Using the definition that the skew is calculated by $ \operatorname{E}\Big[\big(\tfrac{X-\mu}{\sigma}\big)^{\!3}\, \Big]$, and the R formula 
sum((x-mean(x))^3)/(length(x) * sd(x)^3)

I can construct a small, arbitrary distribution to make the skewness low. For example, the distribution 
x = c(1, 3.122, 5, 4, 1.1) 

yields a skew of $-5.64947\cdot10^{-5}$. But this is a small sample and moreover the deviation from symmetry is not large. So, is it possible to construct a larger distribution with one peak that is highly asymmetrical but still has a skewness of nearly zero? 
 A: Consider a distribution on the positive half of the real line which increases linearly from 0 to the mode and then is exponential to the right of the mode, but is continuous at the mode.
This could be called a triangular-exponential distribution (though it does often look a bit like a shark fin).
Let $\theta$ be the location of the mode and $\lambda$ be the rate parameter of the exponential.
As $\lambda\theta$ increases the distribution becomes progressively less skew. As $\lambda\theta$ increases past $\approx 6.15$ the third moment crosses from positive to negative:

Brizzi (2006)$^{[1]}$ refers to this family of distributions as the "two-face" distribution, and discusses this crossover point where the third-moment-skewness is zero. von Hippel (2005)$^{[2]}$ presents an example that's almost at that crossover point here
The thread Non-normal distributions with zero skewness and zero excess kurtosis? has some asymmetric examples, including a small discrete example and another continuous unimodal one:

Discrete unimodal distributions - or equivalently, samples - with zero skewness are quite easy to construct, of large or small size.
Here's an example, which you can treat as a sample or (by dividing the raw frequencies by 3000) as a pmf (the 'x' values are the values taken, the 'n' are the number of times that value occurs in the sample):
x:  -2   -1    0    1    2    3    4    5    6    7    8    9   10
n: 496  498  562 1434    2    1    1    1    1    1    1    1    1


This example is built up from 3-point distributions:
x:          -2              1                  c
n:   c(c-1)(c+1)/6     c(c-1)(c+1)/3 - c       1

across various values of $c$ between 3 and 10. This parameterized (by $c$) 3-point "atom" has $\sum_i n_ix_i =0$ and $\sum_i n_ix_i^3 =0$, which in
turn means that mixtures across various choices of $c$ have zero skewness. (You can't make anything smaller than a distribution across three points that has asymmetry and third central moment zero. A collection of simple pieces over only a few points, such as these make neat building blocks from which larger structures may be made.)
There are all manner of other such "atoms" one can construct, but this example uses only this one kind. To some combination of atoms such as these is added a few symmetrically placed values to fill in remaining holes and guarantee unimodality without destroying the structure of mean and third moment.
$[1]$ Brizzi, M. (2006),
"A Skewed Model Combining Triangular and Exponential Features:
The Two-faced Distribution and its Statistical Properties"
Austrian Journal of Statistics, 35:4, p455–462
http://www.stat.tugraz.at/AJS/ausg064/
$[2]$ von Hippel, P. T. (2005),
"Mean, Median, and Skew: Correcting a Textbook Rule"
Journal of Statistics Education Volume 13, Number 2,
http://ww2.amstat.org/publications/jse/v13n2/vonhippel.html
A: Consider discrete distributions.  One that is supported on $k$ values $x_1, x_2,\ldots, x_k$ is determined by non-negative probabilities $p_1, p_2,\ldots, p_k$ subject to the conditions that (a) they sum to 1 and (b) the skewness coefficient equals 0 (which is equivalent to the third central moment being zero).  That leaves $k-2$ degrees of freedom (in the equation-solving sense, not the statistical one!).  We can hope to find solutions that are unimodal.
To make the search for examples easier, I sought solutions supported on a small symmetrical vector $\mathbf{x}=(-3,-2,-1,0,1,2,3)$ with a unique mode at $0$, zero mean, and zero skewness.  One such solution is $(p_1, \ldots, p_7) = (1396, 3286, 9586, 47386, 8781, 3930, 1235)/75600$.

You can see it is asymmetric.
Here's a more obviously asymmetric solution with $\mathbf{x} = (-3,-1,0,1,2)$ (which is asymmetric) and $p = (1,18, 72, 13, 4)/108$:

Now it's obvious what's going on: because the mean equals $0$, the negative values contribute $(-3)^3=-27$ and $18 \times (-1)^3=-18$ to the third moment while the positive values contribute $4\times 2^3 = 32$ and $13 \times 1^3 = 13$, exactly balancing the negative contributions.  We can take a symmetric distribution about $0$, such as $\mathbf{x}=(-1,0,1)$ with $\mathbf{p}=(1,4,1)/6$, and shift a little mass from $+1$ to $+2$, a little mass from $+1$ down to $-1$, and a slight amount of mass down to $-3$, keeping the mean at $0$ and the skewness at $0$ as well, while creating an asymmetry.  The same approach will work to maintain zero mean and zero skewness of a continuous distribution while making it asymmetric; if we're not too aggressive with the mass shifting, it will remain unimodal.

Edit: Continuous Distributions
Because the issue keeps coming up, let's give an explicit example with continuous distributions.  Peter Flom had a good idea: look at mixtures of normals.  A mixture of two normals won't do: when its skewness vanishes, it will be symmetric.  The next simplest case is a mixture of three normals.
Mixtures of three normals, after an appropriate choice of location and scale, depend on six real parameters and therefore should have more than enough flexibility to produce an asymmetric, zero-skewness solution.  To find some, we need to know how to compute skewnesses of mixtures of normals.  Among these, we will search for any that are unimodal (it is possible there are none).
Now, in general, the $r^\text{th}$ (non-central) moment of a standard normal distribution is zero when $r$ is odd and otherwise equals $2^{r/2}\Gamma\left(\frac{1-r}{2}\right)/\sqrt{\pi}$.  When we rescale that standard normal distribution to have a standard deviation of $\sigma$, the $r^\text{th}$ moment is multiplied by $\sigma^r$.  When we shift any distribution by $\mu$, the new $r^\text{th}$ moment can be expressed in terms of moments up to and including $r$.  The moment of a mixture of distributions (that is, a weighted average of them) is the same weighted average of the individual moments.  Finally, the skewness is zero exactly when the third central moment is zero, and this is readily computed in terms of the first three moments.
This gives us an algebraic attack on the problem.  One solution I found is an equal mixture of three normals with parameters $(\mu, \sigma)$ equal to $(0,1)$, $(1/2,1)$, and $(0, \sqrt{127/18}) \approx (0, 2.65623)$.  Its mean equals $(0 + 1/2 + 0)/3 = 1/6$.  This image shows the pdf in blue and the pdf of the distribution flipped about its mean in red.  That they differ shows they are both asymmetric.  (The mode is approximately $0.0519216$, unequal to the mean of $1/6$.)  They both have zero skewness by construction.

The plots indicate these are unimodal.  (You can check using Calculus to find local maxima.)
A: Here is one I found at https://www.qualitydigest.com/inside/quality-insider-article/problems-skewness-and-kurtosis-part-one.html# which I find nice and reproduced in R: an inverse Burr or Dagum distribution with shape parameters $k=0.0629$ and $c=18.1484$:
$$g(x) = ckx^{-(c+1)}[1+x^{-c}]^{-(k+1)}$$
It has mean 0.5387, standard deviation 0.2907, skewness 0.0000, and kurtosis 2.0000. The source also calls it the "elephant distribution":

My reproduction in R was created with
library(actuar)
library(knotR)

# a nonsymmetric distribution with zero skewness
# see https://www.qualitydigest.com/inside/quality-insider-article/problems-skewness-and-kurtosis-part-one.html#

c <- 18.1484
k <- 0.0629

x <- seq(0,1.5,by=.0001)

elephant.density <- dinvburr(x, k, c)
plot(x,elephant.density, type="l")
polygon(c(min(x),x),c(min(elephant.density),elephant.density), col="grey")
points(0.8,0.8, pch=19, cex=2)

# "ears" created via https://www.desmos.com/calculator/cahqdxeshd
ear.x <- c(0.686, 0.501, 0.42, 0.68)
ear.y <- c(0.698, 0.315, 1.095, 0.983)

myseg(bezier(cbind(ear.x, ear.y)), type="l")

EX <- gamma(k+1/c)*gamma(1-1/c)/gamma(k) # see p6 of https://wwz.unibas.ch/uploads/tx_x4epublication/23_07.pdf
EX2 <- gamma(k+2/c)*gamma(1-2/c)/gamma(k)
EX3 <- gamma(k+3/c)*gamma(1-3/c)/gamma(k)
(skewness <- (EX3 - 3*EX*(EX2-EX^2)-EX^3)/(EX2-EX^2)^(3/2)) # zero to three digits: 0.0003756196

As this output shows, skewness is not quite zero to four digits for these parameter values. Here is a little optimizer for $k$ and $c$:
   # optimize skewness a bit further
    skewval <- 1

while (skewval > 10^(-10)){
  optskew.k <- uniroot(skewness.fun, lower = k*.95, upper = k*1.1, tol=skewval^2, c=c)
  skewval <- optskew.k$f.root
  k <- optskew.k$root

  optskew.c <- uniroot(skewness.fun, lower = c*.95, upper = c*1.1, tol=skewval^2, k=k)
  skewval <- optskew.c$f.root
  c <- optskew.c$root
}

yielding
> print(c)
[1] 18.89306

> print(k)
[1] 0.05975542

> print(skewval)
[1] -1.131464e-15

A: For zero skewness, we need
$$
\operatorname{E}\Big[\big(\tfrac{X-\mu}{\sigma}\big)^{\!3}\, \Big] = 0
$$
or, equivalently,
$$
\operatorname{E}\Big[\big(\tfrac{X-\mu}{\sigma}\big)^{\!3}\, \Big | X \leq \mu \Big] + \operatorname{E}\Big[\big(\tfrac{X-\mu}{\sigma}\big)^{\!3}\, \Big | X \gt \mu \Big] = 0.
$$
Now, for given mean and variance, pick any two distributions $Y$ and $Z$ with zero mass on the right side of $\mu$ and
$$
\operatorname{E}\Big[\big(\tfrac{Y-\mu}{\sigma}\big)^{\!3}\, \Big] = \operatorname{E}\Big[\big(\tfrac{Z-\mu}{\sigma}\big)^{\!3}\, \Big]
$$
and define $X$ to match $Y$ if left of $\mu$ and $(\mu - Z)$ otherwise. (Don't know the exact notation for this, anyone care to help?)
The resulting distribution will be unimodal if the PDFs of $Y$ and $Z$ are increasing at the left of $\mu$ (in addition to being zero at the right of $\mu$).
A: The following discrete distribution is asymmetric and has null skewness: Prob(-4)=1/3, Prob(1)=1/2, Prob(5)=1/6. I found it in the paper of Doric et al., Qual Quant (2009) 43:481-493; DOI 10.1007/s11135-007-9128-9
A: Sure. Try this:
skew= function (x, na.rm = FALSE) 
 {
    if (na.rm)    x <- x[!is.na(x)]             #remove missing values
    sum((x - mean(x))^3)/(length(x) * sd(x)^3)  #calculate skew   
 }

set.seed(12929883) 
x = c(rnorm(100, 1, .1), rnorm(100, 3.122, .1), rnorm(100,5, .1), rnorm(100, 4, .1), rnorm(100,1.1, .1))

 skew(x)
 plot(density(x))

(You did the hard stuff already!)
