Characteristic features of seismicity with long recurrence intervals can be deduced from analysis of historical earthquake records that inherently suffer from uncertainty in the event locations and magnitudes. A novel method to determine the event epicenters and magnitudes jointly from seismic intensities is proposed. The probability for a set of event epicenter and magnitude is assessed by accounting the fitness between the observed and reference seismic intensities, spatial‐occurrence probability based on seismicity density distribution, and temporal‐occurrence probability from the Gutenberg–Richter magnitude–frequency relationship. A set of event epicenter and magnitude yielding the peak probability is chosen. The validity of the method is tested for both synthetic and instrumental seismic‐intensity data, confirming high accuracy. The method is applied effectively to historical events with written seismic damage records. It is found that the errors generally decrease with increasing number and azimuthal coverage of seismic‐intensity data, and increase with epicentral distances. The method appears to be promising for historical earthquakes of which source properties are poorly known. The method is applicable for assessment of the properties of long‐period seismicity, which is crucial for assessment of potential seismic hazards.