A simple geometrical model for statistical mechanical analysis of cooperative binding of monoionic ligands to globular proteins was developed. It is assumed that the ligand has the charged group, which is able to bind the oppositely charged residues on protein surface by electrostatic interactions, and the hydrophobic region, which causes cooperativity by the interactions between bound ligands. The bound Ligand together with the side chain of the ionic residue is assumed to be mobile over the protein surface within the tangential half sphere with the center at the beta-carbon of the side chain and the radius of the effective ligand length. The potential around a bound ligand was approximated by the cylindrical square-well potential with an exclusion core along the center line. The shape of protein was assumed to be spherical. The distance between charged residues on the protein surface was taken from crystallographic data, Grand partition function was calculated as the sum of the terms of all possible bound states. In each term, the electrostatic free energy was taken into account as the Debye-Huckel type electrostatic potential energy for the distribution of isolated charges, and the partition function due to hydrophobic and/or stacking interactions was calculated as the product of pairwise interactions. The model was applied to the binding of an anionic azo dye or anionic surfactants to hen egg white lysozyme at several experimental conditions. The calculated binding isotherms well-reproduced the experimental data. The values of the model parameters estimated were consistent with the ligand size and the magnitude of hydrophobic interaction. Other detailed information such as species fraction, variance, and bound fraction of each site was obtained. The results show the wide applicability of the present model theory.