In this report, I use the pgmpy package to build a bayesian network which name "asia net". I will decide the probability of each disease given some evidence, the detailed structure is shown in the figure below.
'asia' net structure
The CPD values for this model, made of binary variables, are given in the Table below:
the CPD values
For the first step, I use the TabularCPD class to define the CPDs.
Then we need to connect element relationships into the bayesian model.
# Define edges of the asia model
# [Configure element relationships to the model]
asia_model = BayesianModel([('asia', 'tub'),
('smoke', 'lung'),
('smoke', 'bron'),
('tub', 'either'),
('lung', 'either'),
('bron', 'dysp'),
('either', 'xray'),
('either', 'dysp')])
Now let's load the data onto the network.
# Associate the parameters with the model structure.
# [Load the data onto the network]
asia_model.add_cpds(cpd_asia, cpd_smoking, cpd_tuberculosis, cpd_lungCancer, cpd_bronchitis,
cpd_either, cpd_dyspnoea, cpd_Xray)
After setting up the network, we use the variable elimination method to continuously sum and eliminate the variables in the joint probability, and finally obtain the edge distribution.
# Find exact solution if args.exact is True:
# [Perform variable elimination and deduce the predicted value]
if args.exact is True:
print('\n------[variable elimination]------')
asia_infer = VariableElimination(asia_model)
exact = asia_infer.query(
variables=args.qVars,
evidence=evidence)
print(exact)
Let's see the result for exact inference using the variable elimination technique.
Finally I try to sample the model by gibbs sampling which is a Markov Monte Carlo algorithm used in statistics to approximately extract a sequence of samples from a multivariable probability distribution when direct sampling is difficult. The sequence can be used to approximate joint distributions, edge distributions of partial variables, or to compute integrals. Some variables may be known variables, so there is no need to sample them.
# Find approximate solution and cross entropy if args.gibbs is True:
if args.gibbs is True:
gibbs_chain = GibbsSampling(asia_model)
samples = gibbs_chain.sample(evidence=evidence, size=N)
print('\n------[GibbsSampling]------')
print(samples)
Let's see the result by gibbs sampling(the sample size = 1000)