We investigate a network of degenerate optical parametric oscillators (DOPOs) as a model of the coherent Ising machine, an architecture for solving Ising problems. The network represents the interaction in the Ising model, which is a generalization of a previously proposed one for the two-DOPO case. Dynamics of the DOPOs is described by the Fokker-Planck equation in the positive P representation. We obtain approximate distribution of steady states for arbitrary Ising problems under some ansatz. Using the method of statistical mechanics, we analytically demonstrate that the most probable states in a particular range of the parameters correspond to the true optimal states for two rather simple problems, i.e., fully connected ferromagnetic coupling without or with binary random fields. In particular, for the random-field problem, the distribution correctly detects the phase transition that occurs in the target Ising model with varying the magnitude of the fields.