<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Note
Step1:  dtm.csv contains a document-term matrix representation of the words used in Econtalk transcripts. The columns of the matrix correspond to the words in vocab.txt. The rows in the matrix correspond to the show urls in urls.txt.
Our LDA implementation takes input data as a list of lists of hashable objects (typically words). We can use a utility function to convert the document-term matrix to the list of tokenized documents.
Step3:  Here's a utility method to get the title of a webpage that we'll use later.
Step4:  Let's set up our model. First we created a model definition describing the basic structure of our data. Next we initialize an MCMC state object using the model definition, documents, random number generator, and hyper-parameters.
Step5:  When we first create a state object, the words are randomly assigned to topics. Thus, our perplexity (model score) is quite high. After we start to run the MCMC, the score will drop quickly.
Step6:  Run one iteration of the MCMC to make sure everything is working.
Step7:  Now lets run 1000 generations of the MCMC.
Unfortunately, MCMC is slow going.
Step8:  Now that we've run the MCMC, the perplexity has dropped significantly.
Step9:  pyLDAvis projects the topics into two dimensions using techniques described by Carson Sievert.
Step10:  We can extract the term relevance (shown in the right hand side of the visualization) right from our state object. Here are the 10 most relevant words for each topic
Step11:  We could assign titles to each of these topics. For example, Topic 5 appears to be about the foundations of classical liberalism. Topic 6 is obviously Bitcoin and Software.  Topic 0 is the financial system and monetary policy. Topic 4 seems to be generic words used in most episodes; unfortunately, the prevalence of "don" is a result of my preprocessing which splits up the contraction "don't".
We can also get the topic distributions for each document.
Step12:  Topic 5 appears to be about the theory of classical liberalism. Let's find the 20 episodes which have the highest proportion of words from that topic.
Step13:  We could also find the episodes that have notable discussion of both politics AND the financial system.
Step15:  We can look at the topic distributions as projections of the documents into a much lower dimension (16). 
We can try to find shows that are similar by comparing the topic distributions of the documents.
Step16:  Which Econtalk episodes are most similar, in content, to "Mike Munger on the Division of Labor"?
Step17:  How about episodes similar to "Kling on Freddie and Fannie and the Recent History of the U.S. Housing Market"?
Step18:  The model also gives us distributions over words for each topic.
Step19:  We can use this to find the topics a word is most likely to occur in.
Step20:  What topics are most likely to contain the word "Munger" (as in Mike Munger). The number of equal signs indicates the probability the word is generated by the topic. If a topic isn't shown, it's extremely unlikley to generate the word.
Step21:  Where does Munger come up? In discussing the moral foundations of classical liberalism and microeconomics!
How about the word "lovely"? Russ Roberts uses it often when talking about the Theory of Moral Sentiments. It looks like it also comes up when talking about schools.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import pyLDAvis
import json
import sys
import cPickle
from microscopes.common.rng import rng
from microscopes.lda.definition import model_definition
from microscopes.lda.model import initialize
from microscopes.lda import utils
from microscopes.lda import model, runner
from numpy import genfromtxt 
from numpy import linalg
from numpy import array
Explanation: Note: This is best viewed on NBViewer. It is part of a series on Dirichlet Processes and Nonparametric Bayes.
Nonparametric Latent Dirichlet Allocation
Analysis of the topics of Econtalk
In 2003, a groundbreaking statistical model called "Latent Dirichlet Allocation" was presented by David Blei, Andrew Ng, and Michael Jordan.
LDA provides a method for summarizing the topics discussed in a document. LDA defines topics to be discrete probability distrbutions over words. For an introduction to LDA, see Edwin Chen's post.
The original LDA model requires the number of topics in the document to be specfied as a known parameter of the model. In 2005, Yee Whye Teh and others published a "nonparametric" version of this model that doesn't require the number of topics to be specified. This model uses a prior distribution over the topics called a hierarchical Dirichlet process. I wrote an introduction to this HDP-LDA model earlier this year.
For the last six months, I have been developing a Python-based Gibbs sampler for the HDP-LDA model. This is part of a larger library of "robust, validated Bayesian nonparametric models for discovering structure in data" known as Data Microscopes.
This notebook demonstrates the functionality of this implementation.
The Data Microscopes library is available on anaconda.org for Linux and OS X. microscopes-lda can be installed with:
$ conda install -c datamicroscopes -c distributions microscopes-lda
End of explanation
vocab = genfromtxt('./econtalk-data/vocab.txt', delimiter=",", dtype='str').tolist()
dtm = genfromtxt('./econtalk-data/dtm.csv', delimiter=",", dtype=int)
docs = utils.docs_from_document_term_matrix(dtm, vocab=vocab)
urls = [s.strip() for s in open('./econtalk-data/urls.txt').readlines()]
dtm.shape[1] == len(vocab)
dtm.shape[0] == len(urls)
Explanation: dtm.csv contains a document-term matrix representation of the words used in Econtalk transcripts. The columns of the matrix correspond to the words in vocab.txt. The rows in the matrix correspond to the show urls in urls.txt.
Our LDA implementation takes input data as a list of lists of hashable objects (typically words). We can use a utility function to convert the document-term matrix to the list of tokenized documents.
End of explanation
def get_title(url):
    Scrape webpage title
    
    import lxml.html
    t = lxml.html.parse(url)
    return t.find(".//title").text.split("|")[0].strip()
Explanation: Here's a utility method to get the title of a webpage that we'll use later.
End of explanation
N, V = len(docs), len(vocab)
defn = model_definition(N, V)
prng = rng(12345)
state = initialize(defn, docs, prng,
                        vocab_hp=1,
                        dish_hps={"alpha": .6, "gamma": 2})
r = runner.runner(defn, docs, state, )
Explanation: Let's set up our model. First we created a model definition describing the basic structure of our data. Next we initialize an MCMC state object using the model definition, documents, random number generator, and hyper-parameters.
End of explanation
print "randomly initialized model:"
print " number of documents", defn.n
print " vocabulary size", defn.v
print " perplexity:", state.perplexity(), "num topics:", state.ntopics()
Explanation: When we first create a state object, the words are randomly assigned to topics. Thus, our perplexity (model score) is quite high. After we start to run the MCMC, the score will drop quickly.
End of explanation
%%time
r.run(prng, 1)
Explanation: Run one iteration of the MCMC to make sure everything is working.
End of explanation
%%time
r.run(prng, 500)
with open('./econtalk-data/2015-10-07-state.pkl', 'w') as f:
    cPickle.dump(state, f)
%%time
r.run(prng, 500)
with open('./econtalk-data/2015-10-07-state.pkl', 'w') as f:
    cPickle.dump(state, f)
Explanation: Now lets run 1000 generations of the MCMC.
Unfortunately, MCMC is slow going.
End of explanation
print "after 1000 iterations:"
print " perplexity:", state.perplexity(), "num topics:", state.ntopics()
Explanation: Now that we've run the MCMC, the perplexity has dropped significantly.
End of explanation
vis = pyLDAvis.prepare(**state.pyldavis_data())
pyLDAvis.display(vis)
Explanation: pyLDAvis projects the topics into two dimensions using techniques described by Carson Sievert.
End of explanation
relevance = state.term_relevance_by_topic()
for i, topic in enumerate(relevance):
    print "topic", i, ":",
    for term, _ in topic[:10]:
        print term, 
    print 
Explanation: We can extract the term relevance (shown in the right hand side of the visualization) right from our state object. Here are the 10 most relevant words for each topic:
End of explanation
topic_distributions = state.topic_distribution_by_document()
Explanation: We could assign titles to each of these topics. For example, Topic 5 appears to be about the foundations of classical liberalism. Topic 6 is obviously Bitcoin and Software.  Topic 0 is the financial system and monetary policy. Topic 4 seems to be generic words used in most episodes; unfortunately, the prevalence of "don" is a result of my preprocessing which splits up the contraction "don't".
We can also get the topic distributions for each document.
End of explanation
austrian_topic = 5
foundations_episodes = sorted([(dist[austrian_topic], url) for url, dist in zip(urls, topic_distributions)], reverse=True)
for url in [url for _, url in foundations_episodes][:20]:
    print get_title(url), url
Explanation: Topic 5 appears to be about the theory of classical liberalism. Let's find the 20 episodes which have the highest proportion of words from that topic.
End of explanation
topic_a = 0
topic_b = 1
joint_episodes = [url for url, dist in zip(urls, topic_distributions) if dist[0] > 0.18 and dist[1] > 0.18]
for url in joint_episodes:
    print get_title(url), url
Explanation: We could also find the episodes that have notable discussion of both politics AND the financial system.
End of explanation
def find_similar(url, max_distance=0.2):
    Find episodes most similar to input url.
    
    index = urls.index(url)
    for td, url in zip(topic_distributions, urls):
        if linalg.norm(array(topic_distributions[index]) - array(td)) < max_distance:
            print get_title(url), url
Explanation: We can look at the topic distributions as projections of the documents into a much lower dimension (16). 
We can try to find shows that are similar by comparing the topic distributions of the documents.
End of explanation
find_similar('http://www.econtalk.org/archives/2007/04/mike_munger_on.html')
Explanation: Which Econtalk episodes are most similar, in content, to "Mike Munger on the Division of Labor"?
End of explanation
find_similar('http://www.econtalk.org/archives/2008/09/kling_on_freddi.html')
Explanation: How about episodes similar to "Kling on Freddie and Fannie and the Recent History of the U.S. Housing Market"?
End of explanation
word_dists = state.word_distribution_by_topic()
Explanation: The model also gives us distributions over words for each topic.
End of explanation
def bars(x, scale_factor=10000):
    return int(x * scale_factor) * "="
def topics_related_to_word(word, n=10):
    for wd, rel in zip(word_dists, relevance):
        score = wd[word]
        rel_words = ' '.join([w for w, _ in rel][:n]) 
        if bars(score):
            print bars(score), rel_words
Explanation: We can use this to find the topics a word is most likely to occur in.
End of explanation
topics_related_to_word('munger')
Explanation: What topics are most likely to contain the word "Munger" (as in Mike Munger). The number of equal signs indicates the probability the word is generated by the topic. If a topic isn't shown, it's extremely unlikley to generate the word.
End of explanation
topics_related_to_word('lovely')
Explanation: Where does Munger come up? In discussing the moral foundations of classical liberalism and microeconomics!
How about the word "lovely"? Russ Roberts uses it often when talking about the Theory of Moral Sentiments. It looks like it also comes up when talking about schools.
End of explanation
<END_TASK>