An algorithm is developed to model the three-dimensional distribution function of a sample of stars using only measurements of each star's two-dimensional tangential velocity. The algorithm works with missing data: it reconstructs the three-dimensional distribution from data (velocity measurements) that all have one dimension that is unmeasured (the radial direction). It also accounts for covariant measurement uncertainties on the tangential components. The algorithm is applied to tangential velocities measured in a kinematically unbiased sample of 11,865 stars taken from the Hipparcos catalog, chosen to lie on the main sequence and have well-measured parallaxes. The local stellar distribution function of each of a set of 20 color-selected subsamples is modeled as a mixture of two three-dimensional Gaussian ellipsoids of arbitrary relative responsibility. In the fitting, one Gaussian (the halo) is fixed at the known mean and variance tensor of the Galaxy halo, and the other (the disk) is allowed to take an arbitrary mean and an arbitrary variance tensor. The mean and variance tensors (commonly known as the velocity ellipsoid) of the disk distribution are both found to be strong functions of stellar color, with long-lived populations showing larger dispersion, slower mean rotation velocity, and smaller vertex deviation than short-lived populations. The local standard of rest (LSR) is inferred in the usual way, and the Sun's motion relative to the LSR is found to be (U, V, W)☉ = (10.1, 4.0, 6.7) ± (0.5, 0.8, 0.2) km s-1. Artificial data sets are made and analyzed, with the same error properties as the Hipparcos data, to demonstrate that the analysis is unbiased. The results are shown to be insensitive to the assumption that the distributions are Gaussian.