In this paper, we address the identification problem for the systems characterized by linear time-invariant dynamics with bilinear observation models. More precisely, we consider a suitable parametric description of the system and formulate the identification problem as the estimation of the parameters defining the mathematical model of the system using the observed input-output data. To this end, we propose two probabilistic frameworks. The first framework employs the Maximum Likelihood (ML) approach, which accurately finds the optimal parameter estimates by maximizing a likelihood function. Subsequently, we develop a tractable first-order method to solve the optimization problem corresponding to the proposed ML approach. Additionally, to further improve tractability and computational efficiency of the estimation of the parameters, we introduce an alternative framework based on the Expectation--Maximization (EM) approach, which estimates the parameters using an appropriately designed cost function. We show that the EM cost function is invex, which ensures the existence and uniqueness of the optimal solution. Furthermore, we derive the closed-form solution for the optimal parameters and also prove the recursive feasibility of the EM procedure. Through extensive numerical experiments, the practical implementation of the proposed approaches is demonstrated, and their estimation efficacy is verified and compared, highlighting the effectiveness of the methods to accurately estimate the system parameters and their potential for real-world applications in scenarios involving bilinear observation structures.