Learning structure of the bayesian network, and for that there is two main ways: a. By computing an indépendant test between all the variables (with different depth, see Pearl and Company for that) bFor this idea the most simple way is to consider you have a full network (i don't know the english word, but you begin with a network where all nodes are connected, a n-clique). Let's say you have $X_1,X_2,\dots,X_K$ variables. You begin with the test for depth 0, you test if $X_i$ and $X_j$ are independent, for all $i,j<K$. You delete all the edges where you find independence. After you do for depth 1, you test the independence (with chi2 for example) of $X_i$, $X_j$ knowing $X_l$ for all $i,j,l<K$, you need to test only the variables on the neighborhood (for $X_i$ you take the neighborhood of this variable and you test all the triplets inside). After you do for depth 2, 3, 4, 5 and so on. You stop when there is no possible test to do (for example if after depth 0 you find they are all independent, then its finished).
b. By score, you compute the likelihood for example. Then you add or delete or reverse an arc then you compute the likelihood. If the change improve the likelihood then you are happy, if not you forget the change. You loop till no improvement. For the score you can choose many ways (mainly for the penalization).
Often this method begins by creating a tree. For that you compute all the weight between each pair of variable (K(K-1) computations to do, you can do parallelization). Then you order the edge according the weight and you find the minimal tree inside. After you can do what i said, you take a score (for example likelihood penalized with AIC). Its better to take a score which is such that you can compute the new score in a easy way when you just modify your network by adding or deleting or inverting an edge. The steps are easy, you compute the score, you do a modification (add an edge, delete an edge or invert an edge) then you compute the score. It improves, then its good, it doesn't improve then this change stinks. And you continue till you can't improve the score. This is just the main ideas...
Learning the parameters when the structure is known. In this case you just need to compute the frequencies (because you have multinomial, i assume you have a common bayesian network, and with multinomial the linked distribution in bayesian setting is dirichlet. The maximum a posteriori method show that frequency method and bayesian are then very close).
The use of EM algorithm can only be done in a special setup (where you do some assumption about your missing data), at least in a theoretical point of view. But honestly, bayesian network is already heavy in computation, if you add a EM on top of that it will be just so heavy you can"t use it for streaming application.
So i don"t see why you would compute maximum likelihood estimation to get the parameters knowing the structure (which whatever would give you the fact you need to compute the frequencies... I guess you are just looking to a proof of the formulas.)
If i am not clear (reading myself again i feel its not perfect :D), i can try to add some mathematical notation to make it easier.
I am french, so forgive my english. I have a french book about baysian network i like but i suppose there is no translation. You can try to play with bayesian network python library, the one i know is https://pyagrum.readthedocs.io/en/latest/. Read the notebooks should give you some ideas.