Average bit error and symbol energy
The optimal decision regions are the Voronoi Regions. I dont know, if this is what you are looking after.

import numpy as np

points = np.array([(1,1), (1,-1), (-1,1), (-1,-1), (3,3), (3,0), (3,-3), (0,-3), (-3,-3), (-3,0), (-3,3), (0,3), (5,0), (0,5), (-5,0), (0,-5)])
from scipy.spatial import Voronoi, voronoi_plot_2d
vor = Voronoi(points)
voronoi_plot_2d(vor)
plt.axis('equal'); plt.xlim((-7,7)); plt.ylim((-7,7));




I actually, consider analytically calculating the error probability of this constellation a very hard task, since the regions are not very regular.