<jupyter_start><jupyter_text>Gibbs sampling> Programmatic introduction- toc: true - badges: true- comments: true- author: <NAME>- categories: [ML]  Imports<jupyter_code>import networkx as nx
import itertools
from matplotlib import rc
rc("font", family="serif", size=12)
rc("text", usetex=True)
import daft
import random
import requests
import numpy as np
from IPython.display import HTML
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Defining the network<jupyter_code>network={
    "V": ["Letter", "Grade", "Intelligence", "SAT", "Difficulty"],
    "E": [["Intelligence", "Grade"],
        ["Difficulty", "Grade"],
        ["Intelligence", "SAT"],
        ["Grade", "Letter"]],
    "Vdata": {
        "Letter": {
            "ord": 4,
            "numoutcomes": 2,
            "vals": ["weak", "strong"],
            "parents": ["Grade"],
            "children": None,
            "cprob": {
                "['A']": [.1, .9],
                "['B']": [.4, .6],
                "['C']": [.99, .01]
            }
        },

        "SAT": {
            "ord": 3,
            "numoutcomes": 2,
            "vals": ["lowscore", "highscore"],
            "parents": ["Intelligence"],
            "children": None,
            "cprob": {
                "['low']": [.95, .05],
                "['high']": [.2, .8]
            }
        },

        "Grade": {
            "ord": 2,
            "numoutcomes": 3,
            "vals": ["A", "B", "C"],
            "parents": ["Difficulty", "Intelligence"],
            "children": ["Letter"],
            "cprob": {
                "['easy', 'low']": [.3, .4, .3],
                "['easy', 'high']": [.9, .08, .02],
                "['hard', 'low']": [.05, .25, .7],
                "['hard', 'high']": [.5, .3, .2]
            }
        },

        "Intelligence": {
            "ord": 1,
            "numoutcomes": 2,
            "vals": ["low", "high"],
            "parents": None,
            "children": ["SAT", "Grade"],
            "cprob": [.7, .3]
        },

        "Difficulty": {
            "ord": 0,
            "numoutcomes": 2,
            "vals": ["easy", "hard"],
            "parents": None,
            "children": ["Grade"],
            "cprob":  [.6, .4]
        }
    }
}<jupyter_output><empty_output><jupyter_text>Drawing the Bayesian Network<jupyter_code>pgm = daft.PGM([8, 8], origin=[0, 0])

pgm.add_node(daft.Node('Difficulty',r"Difficulty",2,6,aspect=3))   
pgm.add_node(daft.Node('Intelligence',r"Intelligence",5,6,aspect=3))
pgm.add_node(daft.Node('Grade',r"Grade",3,4,aspect=3))    
pgm.add_node(daft.Node('SAT',r"SAT",6,4,aspect=3))    
pgm.add_node(daft.Node('Letter',r"Letter",4,2,aspect=3))   

for node in network['Vdata']:
    parents=network['Vdata'][node]['parents']
    if parents is not None:        
        for parent in parents:
            pgm.add_edge(parent, node) 
pgm.render()<jupyter_output><empty_output><jupyter_text>Finding the Markov blanket of a node<jupyter_code>def find_markov_blanket(node,network):
    '''
    Find the Markov Blanket of the node in the given network
    Markov Blanket is given by:
    1. The parents of the node
    2. The children of the node
    3. The parents of the children of the node
    '''
    
    mb=[]
    #Finding the parents of the node
    parents=network['Vdata'][node]['parents']
    if parents is not None:
        mb.append(parents)
    
    #Finding children of the node
    children=network['Vdata'][node]['children']
    if children is not None:
        mb.append(children)
        
        #Finding parent of each node
        for child in children:
            parents_child=network['Vdata'][child]['parents']
            if parents is not None:
                mb.append(parents)
                
    #Flattening out list of lists
    mb=list(itertools.chain(*mb)) 
    
    #Removing repeated elements
    mb=list(set(mb))
    return mb
find_markov_blanket('Grade',network)<jupyter_output><empty_output><jupyter_text>Gibbs Sampling Procedures  Assigning a random state to a node in the network<jupyter_code>def pick_random(node,network):
    '''
    Assigns a random state to a given node
    N
    '''
    num_outcomes=network['Vdata'][node]["numoutcomes"]
    random_index=random.randint(0,num_outcomes-1)
    return network['Vdata'][node]["vals"][random_index]
pick_random('SAT',network)<jupyter_output><empty_output><jupyter_text>Pick a random non evidence node to the update in the current iteration<jupyter_code>def pick_random_non_evidence_node(non_evidence_nodes):
    return non_evidence_nodes[random.randint(0,len(non_evidence_nodes)-1)]<jupyter_output><empty_output><jupyter_text>Update the value of a node given assignment in previous iteration<jupyter_code>def get_next_value(node, network,simulation):
    parents_current_node_to_update=network['Vdata'][node]['parents']
    if parents_current_node_to_update is None:
        #The node has no parent and we can update it based on the prior
        cumsum=np.cumsum(network['Vdata'][node]["cprob"])    
    else:
        #Find the row corresponding to the values of the parents in the previous iteration
        #NB We need to maintain the order, so we will do it 
        values_parents=[simulation[-1][parent] for parent in parents_current_node_to_update]
        row=network['Vdata'][node]["cprob"][str(values_parents)]
        cumsum=np.cumsum(row)
    choice=random.random()
    index=np.argmax(cumsum>choice)
    return network['Vdata'][node]["vals"][index]<jupyter_output><empty_output><jupyter_text>Main procedure: Iteratively pick up a non evidence node to update<jupyter_code>def gibbs_sampling(network, evidence, niter=2):
    simulation=[]
    nodes=network['V']
    non_evidence_nodes=[node for node in nodes if node not in evidence.keys()]        
    #First iteration random value for all nodes
    d={}
    for node in nodes:
        d[node]=pick_random(node,network)    
    #Put evidence
    for node in evidence:
        d[node]=evidence[node]        
    simulation.append(d.copy())        
    #Now iterate 
    for count in xrange(niter):
        #Pick up a random node to start
        current_node_to_update=pick_random_non_evidence_node(non_evidence_nodes)
        d[current_node_to_update]=get_next_value(current_node_to_update,network,simulation)
        simulation.append(d.copy())
    return simulation<jupyter_output><empty_output><jupyter_text>Illustration 1  Distribution of Letter given that the student is Intelligent<jupyter_code>iterations=int(1e4)
sim=gibbs_sampling(network, {"Intelligence":"high"},iterations)<jupyter_output><empty_output><jupyter_text>Removing first 10% samples<jupyter_code>after_removing_burnt_samples=sim[iterations/10:]
count={val:0 for val in network['Vdata']['Letter']['vals']}<jupyter_output><empty_output><jupyter_text>Finding the distribution of letter<jupyter_code>for assignment in after_removing_burnt_samples:
    count[assignment['Letter']]+=1<jupyter_output><empty_output><jupyter_text>Counts<jupyter_code>count<jupyter_output><empty_output><jupyter_text>Counts to Probabilites<jupyter_code>probabilites={}
for l in count:
    probabilites[l]=count[l]*1.0/(.90*iterations)
probabilites<jupyter_output><empty_output><jupyter_text>Wait a min! What about the marginal distribution of Letter given NO evidence<jupyter_code>iterations=int(1e4)
sim=gibbs_sampling(network, {},iterations)
after_removing_burnt_samples=sim[iterations/10:]
count={val:0 for val in network['Vdata']['Letter']['vals']}
for assignment in after_removing_burnt_samples:
    count[assignment['Letter']]+=1
probabilites_no_evidence={}
for l in count:
    probabilites_no_evidence[l]=count[l]*1.0/(.90*iterations)
probabilites_no_evidence<jupyter_output><empty_output><jupyter_text>How does the evidence about "Intelligent" student affect the quality of letters?<jupyter_code>plt.figure(figsize=(10, 8))
plt.subplot(2,2,1)
plt.bar(range(2),[probabilites['strong'],probabilites['weak']])
plt.xticks([0.5,1.5],['Strong','Weak'])
plt.title('Letter Quality given Intelligent Student')
plt.ylim((0,1.0))
plt.subplot(2,2,2)
plt.bar(range(2),[probabilites_no_evidence['strong'],probabilites_no_evidence['weak']])
plt.xticks([0.5,1.5],['Strong','Weak'])
plt.title('Letter Quality given no prior information')
plt.ylim((0,1.0));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14535022710972986, 'fraction_numerical': 0.014582835285680135, 'mean_word_length': 3.486327077747989, 'pattern_counts': {'":': 48, '<': 62, '<?xml version=': 0, '>': 64, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '376668', 'license': 'Apache-2.0', 'n_tokens_mistral': 2885, 'n_tokens_neox': 2758, 'n_words': 580, 'path': '_notebooks/2014-05-01-gibbs-sampling.ipynb', 'repo_name': 'nipunbatra/blog3', 'size': '8371'}