Understanding the natal kicks, or birth velocities, of neutron stars are essential for understanding the evolution of massive binaries as well as double neutron star formation. We use maximum likelihood methods as published in Verbunt et al. to analyse a new large dataset of parallaxes and proper motions measured by Deller et al. This sample is roughly three times larger than number of measurements available before. For both the complete sample and its younger part (spin-down ages $tau < 3$ Myr), we find that a bimodal Maxwellian distribution describes the measured parallaxes and proper motions better than a single Maxwellian with probability of 99.3 and 95.0 per cent respectively. The bimodal Maxwellian distribution has three parameters: fraction of low-velocity pulsars and distribution parameters $sigma_1$ and $sigma_2$ for low and high-velocity modes. For a complete sample, these parameters are as follows: $42_{-15}^{+17}$ per cent, $sigma_1=128_{-18}^{+22}$ km s$^{-1}$ and $sigma_2 = 298pm 28$ km s$^{-1}$. For younger pulsars, which are assumed to represent the natal kick, these parameters are as follows: $20_{-10}^{+11}$ per cent, $sigma_1=56_{-15}^{+25}$ km s$^{-1}$ and $sigma_2=336pm 45$ km s$^{-1}$. In the young population, $5pm 3$ per cent of pulsars has velocities less than 60 km s$^{-1}$. We perform multiple Monte Carlo tests for the method taking into account realistic observational selection. We find that the method reliably estimates all parameters of the natal kick distribution. Results of the velocity analysis are weakly sensitive to the exact values of scale-lengths of the Galactic pulsar distribution.