On-fault earthquake magnitude distributions are calculated for northern Caribbean faults using estimates of fault slip and regional seismicity parameters. Integer programming, a combinatorial optimization method, is used to determine the optimal spatial arrangement of earthquakes sampled from a truncated Gutenberg-Richter distribution that minimizes the global misfit in slip rates on a complex fault system. Slip rates and their uncertainty on major faults are derived from a previously published GPS block model for the region, with fault traces determined from offshore geophysical mapping and previously published onshore studies. The optimal spatial arrangement of the sampled earthquakes is compared with the 500-year history of earthquake observations. Rupture segmentation of the subduction interface along the Hispaniola-Puerto Rico Trench (PRT) fault and seismic coupling on the PRT fault appear to exert the primary control over this spatial arrangement. Introducing a rupture barrier for the Hispaniola-PRT fault northwest of Mona Passage, based on geophysical and seismicity observations, and assigning a low slip rate of 2 mm/yr on the PRT fault are most consistent with historical earthquakes in the region. The addition of low slip-rate secondary faults as well as segmentation of the Hispaniola and Septentrional strike-slip fault improves the consistency with historical seismicity. An important observation from the modeling is that varying the slip rate on the PRT fault and different segmentation scenarios result in significant changes to the optimal magnitude distribution on faults farther away. In general, optimal on-fault magnitude distributions are more complex and inter-dependent than is typically assumed in probabilistic seismic hazard analysis and probabilistic tsunami hazard analysis.