  We performed Monte Carlo simulations of different properties of pulsar radio
emission, such as: pulsar periods, pulse-widths, inclination angles and rates
of occurrence of interpulse emission (IP). We used recently available large
data sets of the pulsar periods P, the pulse profile widths W and the magnetic
inclination angle alpha. We also compiled the largest ever database of pulsars
with interpulse emission, divided into the double-pole (DP-IP) and the
single-pole (SP-IP) cases. Their distribution on the P - Pdot diagram strongly
suggests a secular alignment of the magnetic axis from the originally random
orientation. We derived possible parent distribution functions of important
pulsar parameters by means of the Kolmogorov-Smirnov significance test using
the available data sets (P, W, alpha and IP), different models of pulsar radio
beam rho = rho(P) as well as different trial distribution functions of pulsar
period and the inclination angles. The best suited parent period distribution
function is the log-normal distribution, although the gamma function
distribution cannot be excluded. The strongest constraint on derived model
distribution functions was the requirement that the numbers of interpulses were
exactly (within 1sigma errors) at the observed level of occurrences. We found
that a suitable model distribution function for the inclination angle is the
complicated trigonometric function which has two local maxima, one near 0 deg
and the other near 90 deg. The former and the latter implies the right rates of
IP occurrence. It is very unlikely that the pulsar beam deviates significantly
from the circular cross-section. We found that the upper limit for the average
beaming factor fb describing a fraction of the full sphere (called also beaming
fraction) covered by a pulsar beam is about 10%. This implies that the number
of the neutron stars in the Galaxy might be underestimated.
