This paper considers the identification of static hysteresis functions which describe phenomena in mechanical systems, piezoelectric actuators and materials. A solution based on a model with a parallel structure of elementary models (with switching) and the Interacting Multiple Model (IMM) approach is proposed. For each elementary model a separate IMM estimator is implemented. The estimated parameters represent a fusion of values from preset grids, weighted by the IMM mode probabilities. The estimated state of each elementary model is a fusion of the estimated states (from the separate Kalman filters) weighted by the IMM probabilities. The nonlinear identification problem is reduced to a linear one. Results from simulation experiments are presented.