Using general identities for difference operators, as well as a technique of symbolic computation and tools from probability theory, we derive very general kth order (k >= 2) convolution identities for Bernoulli and Euler polynomials. This is achieved by use of an elementary result on uniformly distributed random variables. These identities depend on k positive real parameters, and as special cases we obtain numerous known and new identities for these polynomials. In particular we show that the well-known identities of Mild and Matiyasevich for Bernoulli numbers are special cases of the same general formula. (C) 2015 Elsevier Inc. All rights reserved.