Frequency response matrix (FRM) estimation is an important preprocessing step in system identification. This nonparametric step can give a quick insight into the behavior of a dynamic system without making too many assumptions. Sensor failures or faulty communication links make the measurement data go missing. In this paper, a nonparametric method is developed that identifies a multiple-input-multiple-output system from data with missing samples at the noisy outputs. Here, the systems are considered to be excited by arbitrary (random) inputs. The proposed method gives an accurate FRM estimate with its uncertainty and an estimate of the output noise covariance. In addition, the method gives an estimate of the missing time-domain samples and their uncertainties. If the reference signal is known, the proposed method can also handle partially missing data in one or a combination of the following cases: 1) in both inputs and the outputs; 2) in an errors-in-variables framework; and 3) in the presence of a feedback loop.