So, take N=1 and unit masses and G=1. That is, there is a unit mass inside the shell at r < R, using the coordinates in Fig. 2 in the first link that trackstar posted.
Here it is again. There is a single mass located randomly (uniformly) on the shell.
Then, the gravitational field strength at r is a random vector [$]\vec{E}[$] with random magnitude [$]E[$]. Then [$]E[$] has a probability density [$]p(E)[$] with support on the interval [$]\frac{1}{s_{max}^2}[$] to [$]\frac{1}{s_{min}^2}[$] where [$]s_{min} = R -r[$] and [$]s_{max} = R + r[$]. You should be able to work out this density analytically from the Fig. 2 coordinates.
If there is more than one particle, so [$]N \ge 2[$], I would just use a Monte Carlo. Make [$]N[$] draws, distribute the N particles on the surface, calculate the field at the observation point, and repeat for say 10,000 trials to get the pdf and stats.