An efficient quantum Monte Carlo algorithm for the simulation of bosonic systems on a lattice in a grand canonical ensemble is proposed. It is based on the mapping of bosonic models to the spin models in the limit of the infinite total spin quantum number. It is demonstrated how this limit may be taken explicitly in the algorithm, eliminating the systematic errors. The efficiency of the algorithm is examined for the noninteracting lattice boson model and compared with the stochastic series expansion method with the heat-bath-type scattering probability of the random walker.