This paper proposes a new confidence interval for the density of an outcome from a technical system to quantify the uncertainty. Suppose that $n$ random vectors $\{({ X}_i,Y_i)\}_{i=1}^n$ are observed from a technical system, where $X_i\in \Bbb R^d$ is the vector of predictors and $Y_i\in \Bbb R$ is the outcome. The authors are interested in such situations where the sample size $n$ is small because the collection of the $n$ random vector data is rather expensive. The sample size $n=20$ is used in the application. \par The confidence interval proposed is derived by using the random vectors $\{({ X}_i,Y_i)\}_{i=1}^n$ from the technical system, some random vectors $\{({ X}_i,m({ X}_i))\}_{i=n+1}^{n+L}$ from an imperfect simulation model, $Y_i=m({ X}_i)$ for $i=n+1,\dots,n+L$, and other predictor vectors, where $m({ x})$ is a measurable function of ${ x}$. In this derivation, a surrogate estimate $m_n({ x})$ of $m({ x})$ is constructed by fitting a function such as a thin plate spline to $\{({ X}_i,m({ X}_i))\}_{i=n+1}^{n+L}$. It is shown that the goodness-of-fit of the imperfect simulation model affects the upper limit of the confidence interval width. The maximum of absolute errors $\{|Y_i- m_n({ X}_i)|\}_{i=1}^n$ can be used as an evaluation for the goodness-of-fit of the imperfect simulation model. The efficiency of the confidence interval estimation is illustrated in the application to simulated and real data.