Conventional psychometric function (PF) estimation involves fitting a parametric, unidimensional sigmoid to binary subject responses, which is not readily extendible to higher order PFs. This study presents a nonparametric, Bayesian, multidimensional PF estimator that also relies upon traditional binary subject responses. This technique is built upon probabilistic classification (PC), which attempts to ascertain the subdomains corresponding to each subject response as a function of multiple independent variables. Increased uncertainty in the location of class boundaries results in a greater spread in the PF estimate, which is similar to a parametric PF estimate with a lower slope. PC was evaluated on both one-dimensional (1D) and two-dimensional (2D) simulated auditory PFs across a variety of function shapes and sample numbers. In the 1D case, PC demonstrated equivalent performance to conventional maximum likelihood regression for the same number of simulated responses. In the 2D case, where the responses were distributed across two independent variables, PC accuracy closely matched the accuracy of 1D maximum likelihood estimation at discrete values of the second variable. The flexibility and scalability of the PC formulation make this an excellent option for estimating traditional PFs as well as more complex PFs, which have traditionally lacked rigorous estimation procedures.