  We study the problem of estimating linear response statistics under external
perturbations using time series of unperturbed dynamics. Based on the
fluctuation-dissipation theory, this problem is reformulated as an unsupervised
learning task of estimating a density function. We consider a nonparametric
density estimator formulated by the kernel embedding of distributions with
"Mercer-type" kernels, constructed based on the classical orthogonal
polynomials defined on non-compact domains. While the resulting representation
is analogous to Polynomial Chaos Expansion (PCE), the connection to the
reproducing kernel Hilbert space (RKHS) theory allows one to establish the
uniform convergence of the estimator and to systematically address a practical
question of identifying the PCE basis for a consistent estimation. We also
provide practical conditions for the well-posedness of not only the estimator
but also of the underlying response statistics. Finally, we provide a
statistical error bound for the density estimation that accounts for the
Monte-Carlo averaging over non-i.i.d time series and the biases due to a finite
basis truncation. This error bound provides a means to understand the
feasibility as well as limitation of the kernel embedding with Mercer-type
kernels. Numerically, we verify the effectiveness of the estimator on two
stochastic dynamics with known, yet, non-trivial equilibrium densities.
