The Score Function method used to compute first order probabilistic sensitivities is extended in this work to arbitrary-order derivatives included mixed partial derivatives through the use of multicomplex mathematics. Multicomplex mathematics provides an effective and convenient numerical means to compute the high-order kernel functions with respect to natural parameters or moments (mean and standard deviation) obviating the need to analytically determine the kernel functions. Using these numerical kernel functions, high-order derivatives of the response moments or the probability-of-failure with respect to the parameters of the input distributions can be obtained. Numerical results indicate that the high-order probabilistic sensitivities converge with respect to the number of samples at the same rate as standard Monte Carlo estimates. Implementation of multicomplex mathematics is facilitated through the use of the Cauchy-Riemann matrices; therefore, the extension of common engineering probability distributions to matrix form is presented. Nomenclature h perturbation size © 2015. This manuscript version is made available under the Elsevier user license