## 1 Introduction  
Let \(a(x)\) be some unknown continuous multivariate probability density function. Let \(x_{i}=[x_{1},\ldots,x_{in}]\) be the \(i\)th iid sample generated from the distribution. Assume that \(m\) such samples are available; then, a maximum likelihood estimator \(\hat{a}\) will maximize the following expression:  
\[\hat{a}(x)=\prod_{i=1}^{m}\hat{a}(x_{i}) \tag{1}\]  
If \(X\) and \(Y\) are random variables generated from multivariate distribution (\(X\),\(Y\)),  
\[\sum_{i=1}^{m}[y_{i}-\hat{a}(x_{i})]^{2} \tag{2}\]  
where the dimensionality of some response variable \(Y\) is one and an iid sample from the distribution is represented as \((x_{1},y_{1})\),..., \((x_{m},y_{m})\). Then, the least-square non-parametric estimator will minimize the following expression [1]:  
In regression problems, the dependence of \(y\) on \(x\) is considered to be given by some regression function \(\hat{a}(x)\) belonging to some class of functions [2]. In this paper, assuming that some functional model form \(\hat{a}(x)=f(x;\theta_{1},\ldots,\theta_{n})\) is known or a priori established, parameters \(\theta_{1}\),..., \(\theta_{n}\) are computed by minimizing the least-square expression in Equation (2). When the true density \(\alpha(x)\) lies in the parametric class of functions parameterized by the vector \(\theta=[\theta_{1},\ldots,\theta_{m}]\), then, finding parameters by minimizing Equation (2), i.e., maximizing the likelihood, properties such as convergence, consistency, and lack of bias are satisfied [3]. However, when actual density does not lie in the class of parametric functions, there is an acute need for non-parametric estimation [4].  
Methods in Bayesian statistics condition the problem of learning parameters (\(\theta\)) on the dataset \(D\), which is used to learn the parameters. These methods allow for the specification of priors on the parameters as \(p(\theta)\) and require a formal specification of data likelihood function \(p(D|\theta)\), which is conditioned on the parameters \(\theta\). The learning of parameters and confidence bounds on the parameters occurs using Monte Carlo Markov Chain (MCMC) methods to estimate posterior distribution \(p(\theta|D)\) via the following Bayesian rule:  
\[p(\theta|D)\propto p(D|\theta)\times p(\theta) \tag{3}\]  
To reduce any bias in confidence bounds on the parameters, samples from the initial burn-in period are excluded from the computation of confidence bounds.  
The method described in this paper does not require any knowledge of data likelihoods or user-defined priors. The \"population\"-based method directly generates posterior distribution samples \(p(\theta|D)\) so that confidence bounds on parameters can be generated. There is no need to exclude any initial burn-in period samples. More specifically, a genetic algorithm (GA) approach is used to estimate the parameters and their confidence bounds. This method begins with a random sample (i.e., population) estimate \(p^{0}(\theta|D)\) and iterates over certain GA generations \(k\) so that a steady-state distribution of the sample, \(p^{k}(\theta|D)\), minimizing the average least-square error, is obtained. It is shown that from one iteration to the next, the GA forms a Markov Chain, whose transition matrix can be determined using a fitness measure of minimizing the sum of squared error expression (2). Once the population converges, the final population can be used to establish confidence bounds on the regression parameters. The transition matrix is shown to be irreducible and aperiodic, which guarantees both the uniqueness of the steady-state distribution and convergence to a steady-state, regardless of the initial random population \(p^{0}(\theta|D)\).  
The rest of the paper is organized as follows: In Section 2, preliminaries of the GA Markov Chain framework are described. In Section 3, the approach is applied to a dataset, and the results of the experiments are reported. Section 4 concludes the paper with a summary and directions for future work.