How can I find Conditional Probabilties from dataset points of features (random variables)? I am trying to solve the following at work and will dummify for the sake of making it easier to explain myself and getting an answer. My main query is about Step 4 below. But if something is wrong or missing in any other steps, please do correct me.
My system $\Phi$ has $N$ random variables which are actually real-world data. I don't know the Probabilistic relationship between the data. But I did the following to find probabilities. I want to build a Graphical Causal model from my data. Each RV is the node of the graph.
Step 1:
I have binned all the data in some classes for my model. For the sake of simplicity, let's have 2 classes per RV, and $N=4$ i.e. $\{A,B,C,D\}$.
If states of the system are $S_\Phi$, the total number of states that I now have $2^N =2^4 = 16$, such that $S_\Phi \in \{(a^0b^0c^0d^0), (a^0b^0c^0d^1), (a^0b^0c^1d^0), ... (a^1b^1c^1d^1)\}$.
Step 2:
I have calculated individual $P()$, e.g. $P(a^0) = \frac{\text{count all events for }a^0}{\text{total events}} $. This way I made a matrix of Probabilities which I can use later when I need.
Step 3:
Now I find Independent Probability Distribution with the matrix made above $P(A,B,C,D) = P(A).P(B).P(C).P(D)$.  I verified the distribution by getting the sum of all P()s to be $1$ to validate my code.
Step 4:
How can I find the Conditional Probabilities?
I know by formula that , e.g. $P(D=d^0|B=b^1) = \frac{P(B=b^1 \cap D=d^0)}{P(B=b^1)}$. But how do I calculate $P(B=b^1 \cap D=d^0)$. Do I collect all the events when $B=b^1, D=d^0$ and divide it by $16$ (i.e. total number of actual events)? or do I  leave $A, C$ from my total events. In either case the $P(B\cap D)$ is $0.25$. Is it correct? It looks so simple to be true
How can I verify my conditional Probabilities are correct?
Now, I want to make a Probabilistic Graph from these random variables.

First, if you have $N=4$ binary variables $A, B, C, D$, the total number of possible outcomes is not $4^2$ but $2^4$.
Next, in Step 3 you state:
$$
P(A, B, C, D) = P(A)\;P(B)\;P(C)\;P(D),
$$
which would be true if your four random variables were independent. Note, that in this case, your probabilistic graph is just a set of four isolated nodes without any edges between them. If independence is indeed the case, conditioning doesn't change anything, e.g.:
$$
P(D = d^0|B = b^1) = P(D=d^0),
$$
which is the very meaning of "independence": the outcome of $D$ does not depend on the outcome of $B$, so conditioning on $B$ doesn't change anything.
If, however, you cannot presume the above independence, you need to use the formula you stated. E.g.
$$
P(D = d^0|B = b^1) = \frac{P(B=b^1\cap D=d^0)}{P(B=b^1)}.
$$
And you can approximate $P(B=b^1\cap D=d^0)$ by dividing the number of events for which both $B=b^1$ and $D=d^0$ holds, let's denote it by $\#ev(B=b^1\cap D=d^0)$, by the total number of events, let's call it $E$, not by the number of possible outcomes ($2^4$). (E.g. if you toss a coin, you have only two possible outcomes but you can have thousands of events.) I.e.:
$$
\begin{align}
P(D = d^0|B = b^1) &= \frac{P(B=b^1\cap D=d^0)}{P(B=b^1)}\\
    &\approx \frac{\frac{\#ev(B=b^1\cap D=d^0)}{E}}{\frac{\#ev(B=b^1)}{E}}\\
    &= \frac{\#ev(B=b^1\cap D=d^0)}{\#ev(B=b^1)}.\\
\end{align}
$$
Obtaining a probabilistic graph (usually the Bayesian network), is, in general, not trivial. In particular, note that they are usually not unique. You might want to consult libraries that have been developed for this purpose, e.g. bnlearn. But this also only works if you don't have cycles or confounders.
