  We describe the application of a statistical method to estimate submillimeter
galaxy number counts from confusion limited observations by the Balloon-borne
Large Aperture Submillimeter Telescope (BLAST). Our method is based on a
maximum likelihood fit to the pixel histogram, sometimes called 'P(D)', an
approach which has been used before to probe faint counts, the difference being
that here we advocate its use even for sources with relatively high
signal-to-noise ratios. This method has an advantage over standard techniques
of source extraction in providing an unbiased estimate of the counts from the
bright end down to flux densities well below the confusion limit. We
specifically analyse BLAST observations of a roughly 10 sq. deg. map centered
on the Great Observatories Origins Deep Survey South (GOODS-S) field. We
provide estimates of number counts at the three BLAST wavelengths, 250, 350,
and 500 microns; instead of counting sources in flux bins we estimate the
counts at several flux density nodes connected with power-laws. We observe a
generally very steep slope for the counts of about -3.7 at 250 microns and -4.5
at 350 and 500 microns, over the range ~0.02-0.5 Jy, breaking to a shallower
slope below about 0.015 Jy at all three wavelengths. We also describe how to
estimate the uncertainties and correlations in this method so that the results
can be used for model-fitting. This method should be well-suited for analysis
of data from the Herschel satellite.
