Stochastic methods originally devised for geophysical tomography are adapted to the biomedical optical tomography problem. Frequency domain measurements of modulated NIR light are inverted using a Bayesian approximate extended Kalman filter. Minimum variance updates for the linearized problem are calculated from explicit models of the parameters error covariance, the covariance of the system noise, and the measurement error covariance. The method is not iterative per se, but may be applied iteratively to account for strong nonlinearities. Data-driven zonation is used to dynamically reduce the parameterization for improved efficiency, sensitivity, and stability of the inversion. By modeling the parameters as beta distributed random variables, estimates are kept within feasible limits without ad hoc adjustments. In preliminary studies using synthetic domains we have successfully resolved spatially heterogeneous parameters such as absorption, fluorescence lifetime, and quantum efficiency. The method is shown to be much more accurate and computationally efficient than a more traditional Newton-Raphson method. On a 33 by 33 grid, distributed values of a single unknown parameter can be accurately identified in under 2 minutes on a 350 MHz Pentium.