Recently, in order to model systems with uncertainty, probabilistic modeling with conditional probability distribution are widely noticed. Modeling the systems as probability density function enables us to apply strategies considering unsertainty to control systems. One of the most famous modeling methods is Gaussian process regression (GPR). GPR is often used in the field of control systems and its effectiveness is demonstrated in many papers. However, since GPR can represent only Gaussian distribution, GPR cannot always achieve good performances for systems under nonGaussian noise. This paper focuses on kernel density estimation that can approximate arbitrary probability distributions and proposes a probabilistic modeling method for non-Gaussian systems including ones with dynamic characteristics. At that time, a proposed method assumes that measurement is written as a sum of noise and function. Moreover, a prior distribution of the function assumes to be Gaussian. Under the assumption, the proposed method estimates the function as hidden variables by variational Bayes methods. The proposed method can model arbitrary probability density function for single output systems from data. Numerical simulations demonstrate the effectiveness of the proposed method.