text:
# Building a Sentiment Classifier using Scikit-Learn

Acknowledgement: This is derived from https://towardsdatascience.com/building-a-sentiment-classifier-using-scikit-learn-54c8e7c5d2f0.

Image by AbsolutVision @ pixabay.com

>     **Sentiment analysis**, an important area in Natural Language Processing, is the process of automatically detecting affective states of text. Sentiment analysis is widely applied to voice-of-customer materials such as product reviews in online shopping websites like Amazon, movie reviews or social media. It can be just a basic task of classifying the polarity of a text as being positive/negative or it can go beyond polarity, looking at emotional states such as "happy", "angry", etc.

    Here we will build a classifier that is able to distinguish movie reviews as being either positive or negative. For that, we will use [Large Movie Review Dataset v1.0](http://ai.stanford.edu/~amaas/data/sentiment/aclImdb_v1.tar.gz)(2) of IMDB movie reviews.
This dataset contains 50,000 movie reviews divided evenly into 25k train and 25k test. The labels are balanced between the two classes (positive and negative). Reviews with a score <= 4 out of 10 are labeled negative and those with score >= 7 out of 10 are labeled positive. Neutral reviews are not included in the labeled data. This dataset also contains unlabeled reviews for unsupervised learning; we will not use them here. There are no more than 30 reviews for a particular movie because the ratings of the same movie tend to be correlated. All reviews for a given movie are either in train or test set but not in both, in order to avoid test accuracy gain by memorizing movie-specific terms.

## Data preprocessing

    After the dataset has been downloaded and extracted from archive we have to transform it into a more suitable form for feeding it into a machine learning model for training. We will start by combining all review data into 2 pandas Data Frames representing the train and test datasets, and then saving them as csv files: *imdb_train.csv* and *imdb_test.csv*.

    The Data Frames will have the following form:

|text       |label      |
|:---------:|:---------:|
|review1    |0          |
|review2    |1          |
|review3    |1          |
|.......    |...        |
|reviewN    |0          |

    where:
- review1, review2, ... = the actual text of movie review
- 0 = negative review
- 1 = positive review

    But machine learnng algorithms work only with numerical values. We can't just input the text itself into a machine learning model and have it learn from that. We have to, somehow, represent the text by numbers or vectors of numbers. One way of doing this is by using the **Bag-of-words** model(3), in which a piece of text(often called a **document**) is represented by a vector of the counts of words from a vocabulary in that document. This model doesn't take into account grammar rules or word ordering; all it considers is the frequency of words. If we use the counts of each word independently we name this representation a **unigram**. In general, in a **n-gram** we take into account the counts of each combination of n words from the vocabulary that appears in a given document.

    For example, consider these two documents:

d1: "I am learning"

d2: "Machine learning is cool"

The vocabulary of all words encountered in these two sentences is:

v: [ I, am, learning, machine, is, cool ]

    The unigram representations of d1 and d2:

|unigram(d1)|I       |am      |learning|machine |is      |cool    |
|:---------:|:------:|:------:|:------:|:------:|:------:|:------:|
|           |1       |1       |1       |0       |0       |0       |

|unigram(d2)|I       |am      |learning|machine |is      |cool    |
|:---------:|:------:|:------:|:------:|:------:|:------:|:------:|
|           |0       |0       |1       |1       |1       |1       |

    And, the bigrams of d1 and d2 are:

|bigram(d1) |I I     |I am    |I learning|...|machine am|machine learning|...|cool is|cool cool|
|:---------:|:------:|:------:|:--------:|:-:|:--------:|:--------------:|:-:|:-----:|:-------:|
|           |0       |1       |0         |...|0         |0               |...|0      |0        |

|bigram(d2) |I I     |I am    |I learning|...|machine am|machine learning|...|cool is|cool cool|
|:---------:|:------:|:------:|:--------:|:-:|:--------:|:--------------:|:-:|:-----:|:-------:|
|           |0       |0       |0         |...|0         |1               |...|0      |0        |

    Often, we can achieve slightly better results if instead of counts of words we use something called **term frequency times inverse document frequency** (or **tf-idf**). Maybe it sounds complicated, but it is not. Bear with me, I will explain this. The intuition behind this is the following. So, what's the problem of using just the frequency of terms inside a document? Although some terms may have a high frequency inside documents they may not be so relevant for describing a given document in which they appear. That's because those terms may also have a high frequency across the collection of all documents. For example, a collection of movie reviews may have terms specific to movies/cinematography that are present in almost all documents(they have a high **document frequency**). So, when we encounter those terms in a document this doesn't tell much about whether it is a positive or negative review. We need a way of relating **term frequency** (how frequent a term is inside a document) to **document frequency** (how frequent a term is across the whole collection of documents). That is:

$$\begin{align}\frac{\text{term frequency}}{\text{document frequency}} &= \text{term frequency} \cdot \frac{1}{\text{document frequency}} \\ &= \text{term frequency} \cdot \text{inverse document frequency} \\ &= \text{tf} \cdot \text{idf}\end{align}$$

    Now, there are more ways used to describe both term frequency and inverse document frequency. But the most common way is by putting them on a logarithmic scale:

$$tf(t, d) = log(1+f_{t,d})$$
$$idf(t) = log(\frac{1+N}{1+n_t})$$

    where:
$$\begin{align}f_{t,d} &= \text{count of term } \textbf{t} \text{ in document } \textbf{d} \\
N &= \text{total number of documents} \\
n_t &= \text{number of documents that contain term } \textbf{t}\end{align}$$

    We added 1 in the first logarithm to avoid getting $-\infty$ when $f_{t,d}$ is 0. In the second logarithm we added one fake document to avoid division by zero.

    Before we transform our data into vectors of counts or tf-idf values we should remove English **stopwords**(6)(7). Stopwords are words that are very common in a language and are usually removed in the preprocessing stage of natural text-related tasks like sentiment analysis or search.

    Note that we should construct our vocabulary only based on the training set. When we will process the test data in order to make predictions we should use only the vocabulary constructed in the training phase, the rest of the words will be ignored.

    Now, let's create the data frames from the supplied csv files:

```python
import pandas as pd
```

```python
# Read in the training and test datasets from previously created csv files

imdb_train = pd.read_csv('https://raw.githubusercontent.com/RajeevAtla/HTCS-ML-Team-B/master/csv/imdb_train.csv')
imdb_test = pd.read_csv('https://raw.githubusercontent.com/RajeevAtla/HTCS-ML-Team-B/master/csv/imdb_test.csv')
```

```python
# Display information and first few entries from the training and test datasets

pd.set_option('display.max_colwidth', None)

print ("----- Training dataset Info:")
imdb_train.info(verbose=None, buf=None, max_cols=None, memory_usage=None, null_counts=None)
print ("Training dataset Content:")
print(imdb_train.iloc[:5])

print ("\n----- Test dataset Info:")
imdb_test.info(verbose=None, buf=None, max_cols=None, memory_usage=None, null_counts=None)
print ("Test dataset Content:")
print(imdb_test.iloc[:5])

```

    ----- Training dataset Info:
    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 25000 entries, 0 to 24999
    Data columns (total 2 columns):
     #   Column  Non-Null Count  Dtype
    ---  ------  --------------  -----
     0   text    25000 non-null  object
     1   label   25000 non-null  int64
    dtypes: int64(1), object(1)
    memory usage: 390.8+ KB
    Training dataset Content:
                                                                                text  label
    0                                                                                This is what I was expecting when star trek DS9 premiered. Not to slight DS9. That was a wonderful show in it's own right, however it never really gave the fans more of what they wanted. Enterprise is that show. While having a similarity to the original trek it differs enough to be original in it's own ways. It makes the ideas of exploration exciting to us again. And that was one of the primary ingredients that made the original so loved. Another ingredient to success was the relationships that evolved between the crew members. Viewers really cared deeply for the crew. Enterprise has much promise in this area as well. The chemistry between Bakula and Blalock seems very promising. While sexual tension in a show can often become a crutch, I feel the tensions on enterprise can lead to much more and say alot more than is typical. I think when we deal with such grand scale characters of different races or species even, we get some very interesting ideas and television. Also, we should note the performances, Blalock is very convincing as Vulcan T'pol and Bacula really has a whimsy and strength of character that delivers a great performance. The rest of the cast delivered good performances also. My only gripes are as follows. The theme. It's good it's different, but a little to light hearted for my liking. We need something a little more grand. Doesn't have to be orchestral. Maybe something with a little more electronic sound would suffice. And my one other complaint. They sell too many adds. They could fix this by selling less ads, or making all shows two parters. Otherwise we'll end up seeing the shows final act getting wrapped up way too quickly as was one of my complaints of Voyager.      1
    1                                                                                This film reminded me so much of "A History of Violence" which pretended to be a close study of violence and violent behavior but ended up just being nothing short of a cheap action movie masquerading as some thinking film on violence. Dustin Hoffman and his new British bride move to a small English town and encounter endless harassment from the local drunks who do nothing but hang at the pub all day and make trouble. Don't these men have a job? Anyway, Dustin takes all he can take and by the end of the film he holds up in his house and fights off each one of the drunk attackers by such gruesome means as boiling whiskey poured over someone, feet being blown off by a shotgun and someones head getting caught in a bear trap. Funny that someone would have a need for such a large bear trap in a small British town except maybe put a mans head in it. Sam Peckinpah who made the "Wild Bunch" which also covered the topic of blood letting violence in which no one was spared. But it was done with style, and you believed it. Straw Dogs is not believable. First of all the location is wrong and does not work. Why place it in England? I would think maybe in some inner city location or a small town in the American South in the 1930's or something. Second it is not in my view ever really explained clearly why these men are so quick to violence except maybe they got drunk and felt a need to kill Hoffman and rape his wife. Sam Peckinpah missed the mark on this one.      0
    2                                                                                Theodore Rex is possibly the greatest cinematic experience of all time, and is certainly a milestone in human culture, civilization and artistic expression!! In this compelling intellectual masterpiece, Jonathan R. Betuel aligns himself with the great film makers of the 20th century, such as Francis Ford Copola, Martin Scorcese, Orson Welles and Roman Polanski. The special effects are nothing less than breathtaking, and make any work by Spielberg look trite and elementary. At the time of it's release, Theodore Rex was such a revolutionary gem that it raised the bar of film-making to levels never anticipated by film makers. The concept of making not just a motion picture featuring a dinosaur, but adapting an action packed, thrilling detective novel, co-staring a "talking" dinosaur with a post-modern name such as "Theodore", and an existential female police officer changed humanity as we know it. The world could never be the same after experiencing such magnificent beauty. Watching Theodore Rex is much akin to looking into the face of God and hearing Him say "you are my most beloved creation." This is one of the few films that is simply TO DIE FOR!!!      1
    3                                                                                Although I didn't like Stanley & Iris tremendously as a film, I did admire the acting. Jane Fonda and Robert De Niro are great in this movie. I haven't always been a fan of Fonda's work but here she is delicate and strong at the same time. De Niro has the ability to make every role he portrays into acting gold. He gives a great performance in this film and there is a great scene where he has to take his father to a home for elderly people because he can't care for him anymore that will break your heart. I wouldn't really recommend this film as a great cinematic entertainment, but I will say you won't see much bette acting anywhere.      1
    4  'De Grot' is a terrific Dutch thriller, based on the book written by Tim Krabbé. Another of his books, 'Het Gouden Ei' was made into the great Dutch mystery thriller called 'Spoorloos' ('The Vanishing') in 1988. This one is not as good as that thriller (although much better than the American remake also called 'The Vanishing') but there are times it comes close. Especially the opening moments are terrific. We see a man, later we learn his name is Egon Wagter (Fedja van Huêt), coming from a plane in Thailand. When he picks up his bags it is pretty clear that he is smuggling something across the border. These scenes are perfectly directed, photographed and acted. A kind of suspense is created that you would normally not have in an opening scene like this. Later we see how Egon makes his deal in Thailand with a woman, both stating that they have never done anything like this. From this point the movie is constantly flashback and flash-forward. We see how Egon, still as a child (here played by Erik van der Horst), befriends a guy named Axel (as a kid played by Benja Bruijning). We learn how they grew up as friends, sort of, and how Axel (as an adult played by Marcel Hensema) became a criminal. Egon in the meanwhile goes to college and settles with a woman. Around this time he sometimes meets Axel but does not really want anything to do with him. The movie is chronological in a way. It shows Egon and Axel as kids, than as students, young adults, and in their mid-thirties. But from time to time, like I said, the movie goes back to when they were kids and jumps forward again. Every time we see them as kids it explains something that happens when they are adults. Minor spoilers herein. The title means 'The Cave', and it is the cave that gives the movie its happy ending, although it is in fact not that happy. Like the beginning, the ending is terrific. The middle part of the movie is entertaining and in a way it distracts our attention of the first scenes, only to come back at that point in the end. It is the editing that gives the movie its happy ending, although we can say the dramatic ending is happy in a way as well.      1

    ----- Test dataset Info:
    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 25000 entries, 0 to 24999
    Data columns (total 2 columns):
     #   Column  Non-Null Count  Dtype
    ---  ------  --------------  -----
     0   text    25000 non-null  object
     1   label   25000 non-null  int64
    dtypes: int64(1), object(1)
    memory usage: 390.8+ KB
    Test dataset Content:
                                                                                text  label
    0                                                                                Progeny is about a husband and wife who experience time loss while making love. Completely unaware of what this bizarre experience means they try to go on with their lives. The hubby begins questioning the bizarre event and gets help through a very annoying psychiatrist. He comes to believe that aliens are responsible for this lapse in time and that the unborn baby he once thought was his and his wife's actually belongs to the aliens. If ya ask me, this is a great scifi/horror story. Taking a highly questionable real-life scenario involving alien abduction and hybrid breeding is definite thumbs up from this guy. I love all things related to aliens and this story definitely delivered some good ideas. So if you also share an interest in things extraterrestrial, you should be pretty happy with Progeny. At least story-wise anyways. Unfortunately the movie overall is pretty average. With average acting by all actors. Yep, even by the consistently awesome Mr. Dourif, who still does deliver the best performance. Though the black head doctor, delivers his lines really well. There are a few points in the flick where some of the delivery is cringe or laugh worthy, which is fine in my book. I like them cheesy and this had a little bit of some nice stinky cheese, and I mean that in a good way. Anyways, with a less than stellar script you can't really blame all the actors. I especially didn't care for the Mother Hysteria the film went for. She wanted a baby so badly that she'd neglect and dismiss everything her loving husband (who's a doctor!!) said to her. It almost reached a point where you actually didn't care what happened to her. The Progeny is another flick by Brian Yuzna from the icky-sticky film, Society. Again he delivers some slimy effects, and again he delivers a pretty unique tale of horror. If you're into scifi/horror or are a fan of Dourif and or Yuzna films, there's no real reason not to check out this flick if you get the chance. A generous 7 outta 10.      1
    1                                                                                This movie is very similar to Death Warrant with Jean-Claude Van Damme and also has some similarities to Island of Fire with Jackie Chan and I also heard that there is some other very similar action movies, but this film has a much better action than Death Warrant or even Island of Fire (that's right, the Jackie Chan's movie). Rarely American action movies has such a great action sequences, though there was many negative reviews on this film, it easily beats most of the action movies of that time who were more successful. There were many martial art's scenes, David Bradley was fast as Bruce Lee in this film and what else was good, that fighting scenes were much longer than in most of the American martial art's movies. The shoot-out scenes were similar with John Woo's movies, maybe not that good, but still very exiting. There was also many impressive explosions and one great chase scene. I've seen some other David Bradley's movies, but this one, yet is the best in terms of action. OK, this movie has some cheesy moments, but which movie hasn't? The acting was decent, Charles Napier was incredible and his character was real tough. Adam Clark who played Squid and Yuji Okumoto who played the main bad guy were also very good. Other actors acted pretty well too, though the acting isn't important in this type of movies. If you are action movies fan (I mean the real action movies fan, who really can appreciate the good action), than you must see this film.      1
    2                                                                                This is an excellent movie and I would recommend it to everyone. Mr. Drury's acting is top notch as it always is and he blends well with the other actors in the movie. Can't give away any of the suspense or drama found in the movie. Hell to pay is a must see movie!!! The plot was very suspenseful. I would watch this movie over and over again because it has all the elements of a great western movie. It was very authentic in how they displayed the components dealing with this movie which includes the guns, horses, and clothing. The soundtrack is enjoyable and adds flavor to the movie. James Drury has the right touch when picking out a movie to be involved with. This is a another winner for the western genre. !!!!!      1
    3  The 13th and last RKO Falcon film starts with the mutual injunction by Tom Conway as Tom Lawrence alias the Falcon and Ed Brophy as Goldie of "No dames!" whilst they prepare to go on vacation. While you're still wondering what they're going on vacation from as they hadn't had a job since the beginning of the 1st film in 1941 (with Sanders as Gay though and Jenkins as Goldie) they bump into a woman and get dragged into a seedy industrial espionage caper. They promise to help her when her uncle is murdered, by taking an envelope containing the details of a formula to make substitute industrial diamonds to his business colleague in Miami. Suspect everyone here except the cops here who are after Lawrence  and Goldie for the murder. To console himself Goldie keeps paraphrasing travel brochures: "On the coldest day you can always enjoy the warmth of a nice cosy electric chair" for one. Some nice languid atmospheric nightclub scenes rub shoulders with some especially bad behaviour from the baddies. Favourite bit: the dignified game of hide and seek/hunt the thimble the imperturbable and suave Lawrence has with the baddies on the sleeper train. Least favourite bit: the most embarrassing scene in the entire series in the alligator wrestling hut  definitely thrown in for the kids! All in all not the best in the series but yet another entertaining outing, with an overall satisfying plot and many episodes even in this that make me wish they could have gone on for just a few more years as Columbia did with Boston Blackie, although RKO were churning these out faster. Absolutely no sex, not much violence (in fact none at all by today's high standards), and positively no message all make this type of film anathema to serious people who can only regard movies as an art form that must depend on these three pillars. Three Diet Falcon's were made later with John Calvert in the title role, I don't mind them but could never bring myself to count them into the main series, which Tom Conway had made his own by this time. Sad also that it was all downhill after this for Conway, who moved into TV, voice overs and even played Norman Conquest in Park Plaza 605 rather well in 1953. He also developed serious eye and alcohol problems  I don't know if they were linked  wound up poverty stricken and after a spell in hospital in 1967 was found dead in his girlfriend's bed. For us folk that want to at least we still have his 10 entertaining Falcon's plus a number of other worthy, even classic RKO movies from 1942 to 1946 with which to remember him by.      1
    4                                                                                I LOVE Dr WHo SO much! I believe that David Tennant is the best Dr the show has ever had and Billie Piper the Best companion! I liked the way the Dr and Rose had such a connection and a great relationship and the Dr came close a few times to expressing his love for rose! It sadly came to an end after only 2 seasons. I will miss watching rose heaps and think that the show will not be the same without Rose! But David is still there to make me laugh and make me happy to watch him play this fantastic role! I rate this show 110% it is FANTASTIC! The graphics and monsters in this show are wonderful and every storyline is different but somewhat connected and i have actually learned somethings about love, the world and relationships from this show. Therefore it must be one of the most fantastic shows of all time!      1

### Text vectorization

    Fortunately, for the text vectorization part all the hard work is already done in the Scikit-Learn classes `CountVectorizer`(8) and `TfidfTransformer`(5). We will use these classes to transform our csv files into unigram and bigram matrices(using both counts and tf-idf values). (It turns out that if we only use a n-gram for a large n we don't get a good accuracy, we usually use all n-grams up to some n. So, when we say here bigrams we actually refer to uni+bigrams and when we say unigrams it's just unigrams.) Each row in those matrices will represent a document (review) in our dataset, and each column will represent values associated with each word in the vocabulary (in the case of unigrams) or values associated with each combination of maximum 2 words in the vocabulary (bigrams).

    `CountVectorizer` has a parameter `ngram_range` which expects a tuple of size 2 that controls what n-grams to include. After we constructed a `CountVectorizer` object we should call `.fit()` method with the actual text as a parameter, in order for it to learn the required statistics of our collection of documents. Then, by calling `.transform()` method with our collection of documents it returns the matrix for the n-gram range specified. As the class name suggests, this matrix will contain just the counts. To obtain the tf-idf values, the class `TfidfTransformer` should be used. It has the `.fit()` and `.transform()` methods that are used in a similar way with those of `CountVectorizer`, but they take as input the counts matrix obtained in the previous step and `.transform()` will return a matrix with tf-idf values. We should use `.fit()` only on training data and then store these objects. When we want to evaluate the test score or whenever we want to make a prediction we should use these objects to transform the data before feeding it into our classifier.

    Note that the matrices generated for our train or test data will be huge, and if we store them as normal numpy arrays they will not even fit into RAM. But most of the entries in these matrices will be zero. So, these Scikit-Learn classes are using Scipy sparse matrices(9) (`csr_matrix`(10) to be more exactly), which store just the non-zero entries and save a LOT of space.

    We will use a linear classifier with stochastic gradient descent, `sklearn.linear_model.SGDClassifier`(11), as our model. First we will generate and save our data in 4 forms: unigram and bigram matrix (with both counts and tf-idf values for each). Then we will train and evaluate our model for each these 4 data representations using `SGDClassifier` with the default parameters. After that, we choose the data representation which led to the best score and we will tune the hyper-parameters of our model with this data form using cross-validation in order to obtain the best results.

```python
from sklearn.feature_extraction.text import CountVectorizer, TfidfTransformer
```

#### Unigram Counts

```python
# Create a unigram vectorizer and process the training set to generate a list of words.
# Note that unigram processing is set via the ngram_range parameter

unigram_vectorizer = CountVectorizer(ngram_range=(1, 1))
unigram_vectorizer.fit(imdb_train['text'].values)
```

    CountVectorizer()

```python
# Display the length and a few samples of the unigram vectorizer to show the words that have been extracted

print("Number of words found:", len(unigram_vectorizer.get_feature_names()))
print(unigram_vectorizer.get_feature_names()[10000:10100])
print(unigram_vectorizer.get_feature_names()[50000:50100])
```

    Number of words found: 74849
    ['bête', 'bûsu', 'bürgermeister', 'c1', 'c3', 'c3p0', 'c3po', 'c4', 'c57', 'c__p', 'ca', 'caaaaaaaaaaaaaaaaaaaaaaligulaaaaaaaaaaaaaaaaaaaaaaa', 'caalling', 'caan', 'caas', 'cab', 'cabal', 'caballe', 'caballeros', 'caballo', 'cabals', 'cabana', 'cabanne', 'cabaret', 'cabarnet', 'cabbage', 'cabbages', 'cabbie', 'cabby', 'cabel', 'cabell', 'cabells', 'cabin', 'cabinet', 'cabinets', 'cabins', 'cabiria', 'cable', 'cables', 'cabo', 'cabot', 'cabrón', 'cabs', 'caca', 'caccia', 'cache', 'cachet', 'cacho', 'cack', 'cackle', 'cackles', 'cackling', 'caco', 'cacophonist', 'cacophonous', 'cacophony', 'cacoyanis', 'cacoyannis', 'cacti', 'cactus', 'cactuses', 'cad', 'cada', 'cadaver', 'cadaverous', 'cadavers', 'cadavra', 'cadby', 'caddie', 'caddy', 'caddyshack', 'cadena', 'cadence', 'cadences', 'cadet', 'cadets', 'cadfile', 'cadilac', 'cadillac', 'cadillacs', 'cadmus', 'cadre', 'cads', 'caduta', 'cady', 'caesar', 'caesars', 'caeser', 'caetano', 'cafe', 'cafes', 'cafeteria', 'caffari', 'caffeinated', 'caffeine', 'cafferty', 'caffey', 'café', 'cafés', 'cage']
    ['pincher', 'pinchers', 'pinches', 'pinching', 'pinchot', 'pinciotti', 'pine', 'pineal', 'pineapple', 'pineapples', 'pines', 'pinet', 'pinetrees', 'pineyro', 'pinfall', 'pinfold', 'ping', 'pingo', 'pinhead', 'pinheads', 'pinho', 'pining', 'pinjar', 'pink', 'pinkerton', 'pinkett', 'pinkie', 'pinkins', 'pinkish', 'pinko', 'pinks', 'pinku', 'pinkus', 'pinky', 'pinnacle', 'pinnacles', 'pinned', 'pinning', 'pinnings', 'pinnochio', 'pinnocioesque', 'pino', 'pinocchio', 'pinochet', 'pinochets', 'pinoy', 'pinpoint', 'pinpoints', 'pins', 'pinsent', 'pint', 'pinta', 'pinter', 'pintilie', 'pinto', 'pintos', 'pints', 'pinup', 'pioneer', 'pioneered', 'pioneering', 'pioneers', 'piotr', 'pious', 'piovani', 'pip', 'pipe', 'piped', 'pipedream', 'pipeline', 'piper', 'pipers', 'pipes', 'piping', 'pippi', 'pippin', 'pipsqueak', 'piquant', 'piquantly', 'pique', 'piqued', 'piquer', 'piquor', 'piracy', 'pirahna', 'piranha', 'piranhas', 'pirate', 'pirated', 'piraters', 'pirates', 'pirotess', 'pirouette', 'pirouettes', 'pirouetting', 'pirovitch', 'pirro', 'pis', 'pisa', 'pisana']

```python
# Now process the training dataset to get a count of the words extracted earlier

X_train_unigram = unigram_vectorizer.transform(imdb_train['text'].values)
```

```python
# Display the attributes the word count matrix; notice it is huge with 25000 rows since we have 25000 entries
# in the training dataset and 74849 columns since we saw above that we have a vocabulary of 74849 words

print(repr(X_train_unigram))
```

    <25000x74849 sparse matrix of type '<class 'numpy.int64'>'
    	with 3431196 stored elements in Compressed Sparse Row format>

#### Unigram Tf-Idf

```python
# Create a unigram tf-idf vectorizer and load the training set using the word count matrix from earlier

unigram_tf_idf_transformer = TfidfTransformer()
unigram_tf_idf_transformer.fit(X_train_unigram)
```

    TfidfTransformer()

```python
# Now calculate the unigram tf-idf statistics

X_train_unigram_tf_idf = unigram_tf_idf_transformer.transform(X_train_unigram)
```

```python
# Display the attributes the unigram tf-idf matrix; it should be the same size as the unigram matrix above

print(repr(X_train_unigram_tf_idf))
```

    <25000x74849 sparse matrix of type '<class 'numpy.float64'>'
    	with 3431196 stored elements in Compressed Sparse Row format>

#### Bigram Counts

```python
# Create a bigram vectorizer and process the training set to generate a list of bigrams.
# Note that bigram processing is set via the ngram_range parameter and so includes unigrams and bigrams

bigram_vectorizer = CountVectorizer(ngram_range=(1, 2))
bigram_vectorizer.fit(imdb_train['text'].values)
```

    CountVectorizer(ngram_range=(1, 2))

```python
# Display the length and a few samples of the bigram vectorizer to show the bigrams that have been extracted

print("Number of bigrams found:", len(bigram_vectorizer.get_feature_names()))
print(bigram_vectorizer.get_feature_names()[10000:10100])
print(bigram_vectorizer.get_feature_names()[50000:50100])
```

    Number of bigrams found: 1520266
    ['3am but', '3am invesment', '3am it', '3am taped', '3bs', '3bs who', '3d', '3d adventure', '3d although', '3d and', '3d animated', '3d animation', '3d animations', '3d animators', '3d artists', '3d assante', '3d bird', '3d bore', '3d capabilities', '3d cg', '3d cgi', '3d character', '3d companies', '3d computer', '3d dept', '3d disgrace', '3d effect', '3d effects', '3d element', '3d ending', '3d environment', '3d especially', '3d feel', '3d game', '3d games', '3d glasses', '3d grafics', '3d graphics', '3d gravity', '3d had', '3d imax', '3d in', '3d it', '3d just', '3d mario', '3d models', '3d movie', '3d movies', '3d panoramic', '3d plat', '3d probably', '3d programs', '3d sequence', '3d shoot', '3d shooter', '3d shooters', '3d special', '3d splatter', '3d technology', '3d that', '3d trust', '3d was', '3d wave', '3d well', '3d which', '3d with', '3d world', '3d wow', '3dvd', '3dvd collection', '3k', '3k all', '3k attack', '3k but', '3k if', '3k it', '3k not', '3lbs', '3lbs is', '3m', '3m to', '3mins', '3mins jack', '3p', '3p felt', '3p uses', '3p wasn', '3p0', '3p0 r2', '3pm', '3pm and', '3po', '3po and', '3po frank', '3po mind', '3rd', '3rd 4th', '3rd act', '3rd and', '3rd annual']
    ['allen radiation', 'allen randall', 'allen really', 'allen reprises', 'allen reprising', 'allen retains', 'allen retired', 'allen rivkin', 'allen screen', 'allen scripts', 'allen second', 'allen seems', 'allen sees', 'allen set', 'allen she', 'allen should', 'allen showed', 'allen shows', 'allen since', 'allen slowly', 'allen small', 'allen smithee', 'allen soule', 'allen starrer', 'allen stereotypical', 'allen still', 'allen straight', 'allen stuff', 'allen stuttering', 'allen style', 'allen such', 'allen surprisingly', 'allen surrogate', 'allen take', 'allen technique', 'allen tells', 'allen the', 'allen they', 'allen think', 'allen this', 'allen ticks', 'allen to', 'allen treat', 'allen waking', 'allen wannabe', 'allen warning', 'allen was', 'allen we', 'allen when', 'allen whining', 'allen who', 'allen will', 'allen williams', 'allen with', 'allen work', 'allen would', 'allen writes', 'allen writing', 'allen yes', 'allen you', 'allen young', 'allende', 'allende and', 'allende author', 'allende book', 'allende characters', 'allende has', 'allende killed', 'allende magical', 'allende not', 'allende novel', 'allende overthrown', 'allende should', 'allende was', 'allende who', 'allens', 'allens stuff', 'allergic', 'allergic are', 'allergic reaction', 'allergic reactions', 'allergic to', 'allergies', 'allergies kept', 'allergy', 'allergy senator', 'allergy to', 'alles', 'alles action', 'alleviate', 'alleviate her', 'alleviate pain', 'alleviate the', 'alleviate their', 'alleviate visual', 'alley', 'alley 1912', 'alley about', 'alley acting', 'alley again']

```python
# Now generate bigram statistics on the training set

X_train_bigram = bigram_vectorizer.transform(imdb_train['text'].values)
```

```python
# Display the attributes the bigram count matrix; notice it is really huge with 25000 rows since we have 25000 entries
# in the training dataset and 1520266 columns since we saw above that we have 1520266 bigrams

print(repr(X_train_bigram))
```

    <25000x1520266 sparse matrix of type '<class 'numpy.int64'>'
    	with 8689547 stored elements in Compressed Sparse Row format>

#### Bigram Tf-Idf

```python
# Create a bigram tf-idf vectorizer and load the training set using the bigram count matrix from earlier

bigram_tf_idf_transformer = TfidfTransformer()
bigram_tf_idf_transformer.fit(X_train_bigram)
```

    TfidfTransformer()

```python
# Now calculate the bigram tf-idf statistics

X_train_bigram_tf_idf = bigram_tf_idf_transformer.transform(X_train_bigram)
```

```python
# Display the attributes the bigram tf-idf matrix; it should be the same size as the bigram matrix above

print(repr(X_train_bigram_tf_idf))
```

    <25000x1520266 sparse matrix of type '<class 'numpy.float64'>'
    	with 8689547 stored elements in Compressed Sparse Row format>

### Trigram Counts

```python
# Create a trigram vectorizer and process the training set to generate a list of words.
# Note that unigram processing is set via the ngram_range parameter

trigram_vectorizer = CountVectorizer(ngram_range=(1, 3), stop_words={"i", "me", "my", "myself", "we", "our", "ours", "ourselves", "you", "your", "yours", "yourself", "yourselves", "he", "him", "his", "himself", "she", "her", "hers", "herself", "it", "its", "itself", "they", "them", "their", "theirs", "themselves", "what", "which", "who", "whom", "this", "that", "these", "those", "am", "is", "are", "was", "were", "be", "been", "being", "have", "has", "had", "having", "do", "does", "did", "doing", "a", "an", "the", "and", "but", "if", "or", "because", "as", "until", "while", "of", "at", "by", "for", "with", "about", "against", "between", "into", "through", "during", "before", "after", "above", "below", "to", "from", "up", "down", "in", "out", "on", "off", "over", "under", "again", "further", "then", "once", "here", "there", "when", "where", "why", "how", "all", "any", "both", "each", "few", "more", "most", "other", "some", "such", "no", "nor", "not", "only", "own", "same", "so", "than", "too", "very", "s", "t", "can", "will", "just", "don", "should", "now"}) #Generic StopWords
trigram_vectorizer.fit(imdb_train['text'].values)
```

    CountVectorizer(ngram_range=(1, 3),
                    stop_words={'a', 'about', 'above', 'after', 'again', 'against',
                                'all', 'am', 'an', 'and', 'any', 'are', 'as', 'at',
                                'be', 'because', 'been', 'before', 'being', 'below',
                                'between', 'both', 'but', 'by', 'can', 'did', 'do',
                                'does', 'doing', 'don', ...})

```python
# Display the length and a few samples of the trigram vectorizer to show the words that have been extracted

print("Number of words found:", len(unigram_vectorizer.get_feature_names()))
print(trigram_vectorizer.get_feature_names()[10000:10100])
print(trigram_vectorizer.get_feature_names()[50000:50100])
```

    Number of words found: 74849
    ['1920s three cases', '1920s time', '1920s time many', '1920s yet', '1920s yet martians', '1921', '1921 excessively', '1921 excessively affectionate', '1921 fritz', '1921 fritz lang', '1921 irrespective', '1921 irrespective wishes', '1921 mild', '1921 mild jumping', '1921 nicholson', '1921 nicholson tuxedoed', '1921 retired', '1921 retired 1926', '1921 someone', '1921 someone looks', '1921 working', '1921 working steadily', '1922', '1922 account', '1922 account agree', '1922 could', '1922 could see', '1922 dense', '1922 dense fog', '1922 doubt', '1922 doubt best', '1922 dublin', '1922 dublin evoked', '1922 dublin six', '1922 dublin starving', '1922 famous', '1922 famous direct', '1922 film', '1922 film already', '1922 film renoir', '1922 hamdi', '1922 hamdi tell', '1922 hight', '1922 hight black', '1922 learned', '1922 learned lessons', '1922 made', '1922 made toast', '1922 make', '1922 make whatever', '1922 making', '1922 making biggest', '1922 mascot', '1922 mascot 1933', '1922 maybe', '1922 maybe watched', '1922 never', '1922 never seen', '1922 one', '1922 one point', '1922 one thing', '1922 parlor', '1922 parlor prison', '1922 pastor', '1922 pastor open', '1922 screenplay', '1922 screenplay adapted', '1922 silent', '1922 silent film', '1922 still', '1922 still utterly', '1922 testament', '1922 testament dr', '1923', '1923 america', '1923 america time', '1924', '1924 1940', '1924 1940 king', '1924 bauman', '1924 bauman incident', '1924 bring', '1924 bring professional', '1924 death', '1924 death crated', '1924 early', '1924 early 20th', '1924 famous', '1924 famous game', '1924 griffith', '1924 griffith needed', '1924 iron', '1924 iron horse', '1924 kriemhild', '1924 kriemhild margarete', '1924 playing', '1924 playing another', '1924 routing', '1924 routing bonus', '1924 sabu']
    ['accent appears', 'accent appears mish', 'accent ashamed', 'accent ashamed participation', 'accent asian', 'accent asian characters', 'accent aside', 'accent aside serial', 'accent astoundingly', 'accent astoundingly rubbish', 'accent attitude', 'accent attitude pure', 'accent australian', 'accent australian nice', 'accent avoid', 'accent avoid many', 'accent awful', 'accent awful sex', 'accent bacall', 'accent bacall 100', 'accent battle', 'accent battle hulking', 'accent befoul', 'accent befoul film', 'accent beginning', 'accent beginning american', 'accent believe', 'accent believe french', 'accent belting', 'accent belting adelaide', 'accent ben', 'accent ben chaplin', 'accent besides', 'accent besides seems', 'accent better', 'accent better script', 'accent beyond', 'accent beyond insulting', 'accent bizarre', 'accent bizarre although', 'accent bizarre stan', 'accent bland', 'accent bland often', 'accent body', 'accent body language', 'accent bought', 'accent bought dvd', 'accent brad', 'accent brad johnson', 'accent brad pitt', 'accent breaking', 'accent breaking much', 'accent british', 'accent british maybe', 'accent butler', 'accent butler hovis', 'accent buy', 'accent buy passably', 'accent california', 'accent california looks', 'accent canadian', 'accent canadian mountie', 'accent canned', 'accent canned laugh', 'accent cast', 'accent cast assembled', 'accent cause', 'accent cause helped', 'accent certainly', 'accent certainly makes', 'accent character', 'accent character supposed', 'accent chemistry', 'accent chemistry two', 'accent chose', 'accent chose typical', 'accent cinematography', 'accent cinematography editing', 'accent claire', 'accent claire danes', 'accent classic', 'accent classic character', 'accent clearly', 'accent clearly ball', 'accent close', 'accent close major', 'accent clumsy', 'accent clumsy barely', 'accent comes', 'accent comes goes', 'accent compare', 'accent compare friday', 'accent concerned', 'accent concerned must', 'accent conjured', 'accent conjured picture', 'accent control', 'accent control body', 'accent conveniently', 'accent conveniently volunteered']

```python
# Now process the training dataset to get a count of the words extracted earlier

X_train_trigram = trigram_vectorizer.transform(imdb_train['text'].values)
```

### Trigram Tf-Idf

```python
trigram_tf_idf_transformer = TfidfTransformer()
trigram_tf_idf_transformer.fit(X_train_trigram)
```

    TfidfTransformer()

```python
# Now calculate the trigram tf-idf statistics

X_train_trigram_tf_idf = trigram_tf_idf_transformer.transform(X_train_trigram)
```

```python
# Display the attributes the trigram tf-idf matrix; it should be the same size as the trigram matrix above

print(repr(X_train_trigram_tf_idf))
```

    <25000x4696483 sparse matrix of type '<class 'numpy.float64'>'
    	with 8435367 stored elements in Compressed Sparse Row format>

### Try the four different data formats (unigram, bigram with and without tf_idf) on the training set and pick the best

    Now, for each data form we split it into train & validation sets, train a `SGDClassifier` and output the score.

```python
from sklearn.linear_model import SGDClassifier
from sklearn.model_selection import train_test_split
from scipy.sparse import csr_matrix
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from typing import Tuple
```

```python
# Helper function to display confusion matrix

def display_confusion_matrix(y_true, y_pred) -> None:
    cf_matrix = confusion_matrix(y_true, y_pred)
    group_names = ['True Neg','False Pos','False Neg','True Pos']
    group_counts = ["{0:0.0f}".format(value) for value in cf_matrix.flatten()]
    group_percentages = ["{0:.2%}".format(value) for value in cf_matrix.flatten()/np.sum(cf_matrix)]
    labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in zip(group_names,group_counts,group_percentages)]
    labels = np.asarray(labels).reshape(2,2)
    sns.heatmap(cf_matrix, annot=labels, fmt='', cmap='coolwarm')
    plt.show()

```

```python
def train_and_show_scores(X: csr_matrix, y: np.array, title: str) -> Tuple[float, float]:
    X_train, X_valid, y_train, y_valid = train_test_split(
        X, y, train_size=0.75, stratify=y
    )

    clf = SGDClassifier()
    clf.fit(X_train, y_train)
    train_score = clf.score(X_train, y_train)
    valid_score = clf.score(X_valid, y_valid)
    print(f'{title}\nTrain score: {round(train_score, 2)} ; Validation score: {round(valid_score, 2)}')

    train_pred = clf.predict(X_train)
    valid_pred = clf.predict(X_valid)
    print(f'Train precision: {round(precision_score(y_train, train_pred), 2)} ; Validation precision: {round(precision_score(y_valid, valid_pred), 2)}')
    print(f'Train recall: {round(recall_score(y_train, train_pred), 2)} ; Validation recall: {round(recall_score(y_valid, valid_pred), 2)}')
    print(f'Train F1: {round(f1_score(y_train, train_pred), 2)} ; Validation F1: {round(f1_score(y_valid, valid_pred), 2)}')
    print("Train Confusion Matrix: ")
    print(confusion_matrix(y_train, train_pred))
    display_confusion_matrix(y_train, train_pred)
    print("Validation Confusion Matrix: ")
    print(confusion_matrix(y_valid, valid_pred))
    display_confusion_matrix(y_valid, valid_pred)
    print("\n")

    return train_score, valid_score, clf
```

```python
y_train = imdb_train['label'].values
```

```python
uc_train_score, uc_valid_score, uc_clf = train_and_show_scores(X_train_unigram, y_train, '----- Unigram Counts -----')
utfidf_train_score, utfidf_valid_score, utfidf_clf = train_and_show_scores(X_train_unigram_tf_idf, y_train, '----- Unigram Tf-Idf -----')
bc_train_score, bc_valid_score, bc_clf = train_and_show_scores(X_train_bigram, y_train, '----- Bigram Counts -----')
btfidf_train_score, btfidf_valid_score, btfidf_clf = train_and_show_scores(X_train_bigram_tf_idf, y_train, '----- Bigram Tf-Idf -----')
tc_train_score, tc_valid_score, tc_clf = train_and_show_scores(X_train_trigram, y_train, '----- Trigram Counts -----')
ttfidf_train_score, ttfidf_valid_score, ttfidf_clf = train_and_show_scores(X_train_trigram_tf_idf, y_train, '----- Trigram Tf-Idf -----')
```

```python
# Display the previously derived scores for the four scenarios

sns.set_style("whitegrid", {'grid.linestyle': '--'})

print ("Training score for the four approaches:")
ax1 = sns.barplot(
    x= ['Unigram Count', 'Unigram tf-idf', 'Bigram Count', 'Bigram tf-idf'],
    y= [uc_train_score, utfidf_train_score, bc_train_score, btfidf_train_score])
ax1.set(ylim=(0.8, 1.0))
plt.show()

print ("Validation score for the four approaches:")
ax2 = sns.barplot(
    x= ['Unigram Count', 'Unigram tf-idf', 'Bigram Count', 'Bigram tf-idf'],
    y= [uc_valid_score, utfidf_valid_score, bc_valid_score, btfidf_valid_score])
ax2.set(ylim=(0.8, 1.0))
plt.show()
```

    The best data form seems to be **bigram with tf-idf** as it gets the highest validation accuracy: **0.9**; so we will choose it as our preferred approach.

### Testing model

```python
# Transform the test data set into the bigram tf-idf format

X_test = bigram_vectorizer.transform(imdb_test['text'].values)
X_test = bigram_tf_idf_transformer.transform(X_test)
y_test = imdb_test['label'].values
```

```python
# Now evaluate the test data using the previously trained bigram tf-idf classifier

clf = btfidf_clf
score = clf.score(X_test, y_test)

print(f'Score: {round(score, 4)}')

test_pred = clf.predict(X_test)
print(f'Test precision: {round(precision_score(y_test, test_pred), 4)}')
print(f'Test recall: {round(recall_score(y_test, test_pred), 4)}')
print(f'Test F1: {round(f1_score(y_test, test_pred), 4)}')
print("Test Confusion Matrix: ")
print(confusion_matrix(y_test, test_pred))
display_confusion_matrix(y_test, test_pred)
print("\n")
```

```python
#import all necessary packages
from sklearn.linear_model import SGDClassifier
from sklearn.model_selection import train_test_split
from scipy.sparse import csr_matrix
from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score, accuracy_score
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from typing import Tuple
import pandas as pd
from sklearn.feature_extraction.text import CountVectorizer, TfidfTransformer
import sys
from PIL import Image
```

```python
# Read in the training and test datasets from previously created csv files

imdb_train = pd.read_csv('https://raw.githubusercontent.com/RajeevAtla/HTCS-ML-Team-B/master/csv/imdb_train.csv')
imdb_test = pd.read_csv('https://raw.githubusercontent.com/RajeevAtla/HTCS-ML-Team-B/master/csv/imdb_test.csv')
```

```python
# Define the stop words

STOP_WORDS = {"i", "me", "my", "myself", "we", "our", "ours", "ourselves", "you",
              "your", "yours", "yourself", "yourselves", "he", "him", "his", "himself",
              "she", "her", "hers", "herself", "it", "its", "itself", "they", "them",
              "their", "theirs", "themselves", "what", "which", "who", "whom", "this",
              "that", "these", "those", "am", "is", "are", "was", "were", "be", "been",
              "being", "have", "has", "had", "having", "do", "does", "did", "doing",
              "a", "an", "the", "and", "but", "if", "or", "because", "as", "until",
              "while", "of", "at", "by", "for", "with", "about", "against", "between",
              "into", "through", "during", "before", "after", "above", "below", "to",
              "from", "up", "down", "in", "out", "on", "off", "over", "under", "again",
              "further", "then", "once", "here", "there", "when", "where", "why", "how",
              "all", "any", "both", "each", "few", "more", "most", "other", "some", "such",
              "no", "nor", "not", "only", "own", "same", "so", "than", "too", "very",
              "s", "t", "can", "will", "just", "don", "should", "now"}
```

```python
y_train = imdb_train['label'].values

y_test = imdb_test['label'].values

#3 parameters to be changed are loss, alpha, and eta0

#loss is the loss function to use
#possible options are ‘hinge’, ‘log’, ‘modified_huber’, ‘squared_hinge’, ‘perceptron’, ‘squared_loss’, ‘huber’, ‘epsilon_insensitive’, or ‘squared_epsilon_insensitive’
#Google/Wikipedia are your friends here

#alpha is regularization coefficient

#eta0 is initial learning rate

def train_and_show_n_gram (n = 1, bool_stop_words = False, tf_idf = False, loss = 'hinge', alpha = 0.0001, eta0 = 0.0, title = 'unnamed.jpg'):

    #n-gram counts
    if bool_stop_words:
        vectorizer = CountVectorizer(ngram_range = (1, n), stop_words = STOP_WORDS) #only run this if stop words is passed
    else:
        vectorizer = CountVectorizer(ngram_range = (1, n))

    #fit to data
    vectorizer.fit(imdb_train['text'].values)

    X_train_n_gram = vectorizer.transform(imdb_train['text'].values)

    #tf-idf
    if tf_idf:
        tf_idf_transformer = TfidfTransformer()
        tf_idf_transformer.fit(X_train_n_gram)
        train_data = tf_idf_transformer.transform(X_train_n_gram)
    else:
        train_data = X_train_n_gram

    #test-validation split
    X_train, X_valid, Y_train, y_valid = train_test_split(train_data, y_train, train_size = 0.75, stratify = y_train)

    #SGDClassifier with other parameters
    clf = SGDClassifier(loss = loss, alpha = alpha, eta0 = eta0)
    clf.fit(X_train, Y_train)

    #get training and validation scores
    train_score = clf.score(X_train, Y_train)
    valid_score = clf.score(X_valid, y_valid)

    #Evaluate model on test data
    X_test_n_gram = vectorizer.transform(imdb_test['text'].values)

    if tf_idf:
        test_tf_idf_transformer = TfidfTransformer()
        test_tf_idf_transformer.fit(X_test_n_gram)
        test_data = test_tf_idf_transformer.transform(X_test_n_gram)
    else:
        test_data = X_test_n_gram

    y_pred = clf.predict(test_data)

    cf = confusion_matrix(y_test, y_pred)

    group_names = ['True Negatives','False Positives','False Negatives','True Positives']
    group_counts = ["{0:0.0f}".format(value) for value in cf.flatten()]
    group_percentages = ["{0:.2%}".format(value) for value in cf.flatten()/np.sum(cf)]

    labels = [f"{v1}\n{v2}\n{v3}" for v1, v2, v3 in zip(group_names, group_counts, group_percentages)]
    labels = np.asarray(labels).reshape(2, 2)
    sns.heatmap(cf, annot = labels, fmt = '', cmap = 'coolwarm')
    plt.title(label = title) #to identify/differentiate between the confusion matrices

    f1 = f1_score(y_test, y_pred)
    accuracy = accuracy_score(y_test, y_pred)

    #display the test accuracy for each model
    plt.xlabel(f'Test Accuracy: {accuracy}')

    #save and show
    plt.savefig(title)
    plt.show()

    return train_score, valid_score, clf

```

```python
uc_train_score, uc_valid_score, uc_clf = train_and_show_n_gram(n = 1, title = '1.jpg')
utfidf_train_score, utfidf_valid_score, utfidf_clf = train_and_show_n_gram(n = 1, tf_idf = True, title = '1-tf-idf.jpg')
uc_stop_words_train_score, uc_stop_words_valid_score, uc_stop_words_clf = train_and_show_n_gram(n = 1, bool_stop_words = True, title = '1-SW.jpg')
utfidf_stop_words_train_score, utfidf_stop_words_valid_score, utfidf_stop_words_clf = train_and_show_n_gram(n = 1, tf_idf = True, bool_stop_words = True, title = '1-SW-tf-idf.jpg')

bc_train_score, bc_valid_score, bc_clf = train_and_show_n_gram(n = 2, title = '2.jpg')
btfidf_train_score, btfidf_valid_score, btfidf_clf = train_and_show_n_gram(n = 2, tf_idf = True, title = '2-tf-idf.jpg')
bc_stop_words_train_score, bc_stop_words_valid_score, bc_stop_words_clf = train_and_show_n_gram(n = 2, bool_stop_words = True, title = '2-SW.jpg')
btfidf_stop_words_train_score, btfidf_stop_words_valid_score, btfidf_stop_words_clf = train_and_show_n_gram(n = 2, tf_idf = True, bool_stop_words = True, title = '2-SW-tf-idf.jpg')

tc_train_score, tc_valid_score, tc_clf = train_and_show_n_gram(n = 3, title = '3.jpg')
ttfidf_train_score, ttfidf_valid_score, ttfidf_clf = train_and_show_n_gram(n = 3, tf_idf = True, title = '3-tf-idf.jpg')
tc_stop_words_train_score, tc_stop_words_valid_score, tc_stop_words_clf = train_and_show_n_gram(n = 3, bool_stop_words = True, title = '3-SW.jpg')
ttfidf_stop_words_train_score, ttfidf_stop_words_valid_score, ttfidf_stop_words_clf = train_and_show_n_gram(n = 3, tf_idf = True, bool_stop_words = True, title = '3-SW-tf-idf.jpg')

#print(uc_train_score, uc_valid_score)
#print(utfidf_train_score, utfidf_valid_score)

#print(uc_stop_words_train_score, uc_stop_words_valid_score)
#print(utfidf_stop_words_train_score, utfidf_stop_words_valid_score)
```

```python
#the following cells will save the confusion matrices after bunching them together, based on the value of n
images = [Image.open(x) for x in ['1.jpg', '1-tf-idf.jpg', '1-SW.jpg', '1-SW-tf-idf.jpg']]
widths, heights = zip(*(i.size for i in images))

total_width = sum(widths)
max_height = max(heights)

new_im = Image.new('RGB', (total_width, max_height))

x_offset = 0
for im in images:
  new_im.paste(im, (x_offset,0))
  x_offset += im.size[0]

new_im.save('1-final.jpg')
```

```python
images = [Image.open(x) for x in ['2.jpg', '2-tf-idf.jpg', '2-SW.jpg', '2-SW-tf-idf.jpg']]
widths, heights = zip(*(i.size for i in images))

total_width = sum(widths)
max_height = max(heights)

new_im = Image.new('RGB', (total_width, max_height))

x_offset = 0
for im in images:
  new_im.paste(im, (x_offset,0))
  x_offset += im.size[0]

new_im.save('2-final.jpg')
```

```python
images = [Image.open(x) for x in ['3.jpg', '3-tf-idf.jpg', '3-SW.jpg', '3-SW-tf-idf.jpg']]
widths, heights = zip(*(i.size for i in images))

total_width = sum(widths)
max_height = max(heights)

new_im = Image.new('RGB', (total_width, max_height))

x_offset = 0
for im in images:
  new_im.paste(im, (x_offset,0))
  x_offset += im.size[0]

new_im.save('3-final.jpg')
```

```python
#display the training scores across all models
sns.reset_defaults()
sns.plotting_context(font_scale = 1)

plt.figure(figsize=(12, 10))
x_labels = ['1', '1-SW', '1-tf-idf', '1-SW-tf-idf', '2', '2-tf-idf', '2-SW', '2-SW-tf-idf', '3', '3-tf-idf', '3-SW', '3-SW-tf-idf']
values = [uc_train_score, uc_stop_words_train_score, utfidf_train_score, utfidf_stop_words_train_score, bc_train_score, btfidf_train_score, bc_stop_words_train_score, btfidf_stop_words_train_score, tc_train_score, ttfidf_train_score, tc_stop_words_train_score, ttfidf_stop_words_train_score]

ax1 = sns.barplot(
    x = x_labels,
    y = values,
    saturation = 1)

ax1.set(ylim = (0.90, 1.005), title = 'Training Scores', xlabel = 'Model', ylabel = 'Training Score')
plt.setp(ax1, yticks = np.arange(0.9,1.005, 0.0025))
plt.savefig("Training_Scores.jpg", quality = 100)
plt.show()

```

```python
#display the validation scores across all models

sns.reset_defaults()
sns.plotting_context(font_scale = 1)

plt.figure(figsize = (12, 10))

x_labels_test = ['1', '1-SW', '1-tf-idf', '1-SW-tf-idf', '2', '2-tf-idf', '2-SW', '2-SW-tf-idf', '3', '3-tf-idf', '3-SW', '3-SW-tf-idf']
values_test = [uc_valid_score, uc_stop_words_valid_score, utfidf_valid_score, utfidf_stop_words_valid_score, bc_valid_score, btfidf_valid_score, bc_stop_words_valid_score, btfidf_stop_words_valid_score, tc_valid_score, ttfidf_valid_score, tc_stop_words_valid_score, ttfidf_stop_words_valid_score]

ax1 = sns.barplot(
    x = x_labels_test,
    y = values_test,
    saturation = 1)

ax1.set(ylim = (0.7, 0.92), title = 'Validation Scores', xlabel = 'Model', ylabel = 'Validation Score')

plt.setp(ax1, yticks = np.arange(0.7,0.92, 0.005))
plt.savefig("Validation_Scores.jpg", quality = 100)
plt.show()
```

    And we got almost 90% test accuracy. That's not bad for our simple linear model. There are more advanced methods that give better results. The current state-of-the-art on this dataset is **97.42%** (13)

## References

(1)  [Sentiment Analysis - Wikipedia](https://en.wikipedia.org/wiki/Sentiment_analysis)
(2)  [Learning Word Vectors for Sentiment Analysis](http://ai.stanford.edu/~amaas/papers/wvSent_acl2011.pdf)
(3)  [Bag-of-words model - Wikipedia](https://en.wikipedia.org/wiki/Bag-of-words_model)
(4)  [Tf-idf - Wikipedia](https://en.wikipedia.org/wiki/Tf%E2%80%93idf)
(5)  [TfidfTransformer - Scikit-learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.text.TfidfTransformer.html)
(6)  [Stop words - Wikipedia](https://en.wikipedia.org/wiki/Stop_words)
(7)  [A list of English stopwords](https://gist.github.com/sebleier/554280)
(8)  [CountVectorizer - Scikit-learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.text.CountVectorizer.html)
(9)  [Scipy sparse matrices](https://docs.scipy.org/doc/scipy/reference/sparse.html)
(10) [Compressed Sparse Row matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix)
(11) [SGDClassifier - Scikit-learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDClassifier.html)
(12) [RandomizedSearchCV - Scikit-learn documentation](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RandomizedSearchCV.html)
(13) [Sentiment Classification using Document Embeddings trained with
Cosine Similarity](https://www.aclweb.org/anthology/P19-2057.pdf)

max_stars_repo_name: RajeevAtla/HTCS-ML

max_issues_repo_name: RajeevAtla/HTCS-ML

max_forks_repo_name: RajeevAtla/HTCS-ML


text:
# # Exercises for Machine Learning with Python, Lecture 2: *Kernel Ridge Regression and Molecules*
By Dr. Anders Christensen `anders.christensen @ unibas.ch`

---
# ***Important:***

***The answers to the questions must be entered in the online system.***

https://www.chemie1.unibas.ch/~pythonprakt/

*    You should have received the username and password by email.
*    Use you @unibas.ch or @stud.unibas.ch email to login after you enter the username and password for the course.

## Before you begin the exercises:

Again, we begin by downloading a dataset. This dataset is taken from the "QM7" dataset introduced in the lecture.

As there are 7101 molecules in the datasets, all the data is stored in binary format to keep the files small enough for Google Colab. The data format is Numpy's ".npy" format.

The data consists of representations (aka. features) for each of the 7101 molecules, and also the atomization energy of each molecule (i.e. how much energy it takes to break the molecule into its free atoms).

There are two sets of representations, the Coulomb Matrix representation and the Bag-of-Bonds representation, and later you will compare which of the two is best for predicting the atomization energy with Kernel Ridge Regression.

**To-do:**

*    Download the binary datafiles with the "wget" command:

```
# Download features maps and labels:
!wget -O cmat.npy https://www.dropbox.com/s/4eqafm0yr2ywxps/cmat.npy
!wget -O bob.npy https://www.dropbox.com/s/vyiwza2uy4jkczg/bob.npy
!wget -O hof.npy https://www.dropbox.com/s/zy717f8mwxaegff/hof.npy
```

Next, Numpy can directly read .npy files and convert these into Numpy arrays using the `np.load()` function.

**To-do:**

*   Run the code below to read in the features as numpy arrays and print their sizes and types:

```
import numpy as np

# Save 'Coulomb matrices' features in 'cmat_features'
cmat_features = np.load("cmat.npy")

# Save 'Bag-of-Bonds' features in 'bob_features'
bob_features = np.load("bob.npy")

# Save Atomization Energies in 'hof'
atomization_energy = np.load("hof.npy")

# Print the dimensions and type of each numpy array, just to see what you got
print(cmat_features.shape)
print(type(cmat_features))

print(bob_features.shape)
print(type(bob_features))

print(atomization_energy.shape)
print(type(atomization_energy))
```

As you can see above, each Numpy array contains 7101 rows. Each row in the Coulomb-matrix representation has 276 features and each row in the Bag-of-bonds representation has 465 features.

The 7101 heats-of-formation are in units of kcal/mol.

Additionally, you can print all the arrays, just to see what they contain (a lot of numbers).

```
print(cmat_features)
print(bob_features)
print(atomization_energy)
```

# Exercises 2.1: Exploratory Data Analysis

As we've done during the lecture, lets's start by having a look at your data.

### Question 2.1.1: (6.1)

Since we've already looked into the representation during the lecture, we will skip looking at plots for them those for now.

However, to make sure that we trust our data, we still want to plot a histogram of the labels.

**To-do:**
*    Using matplotlib/pyplot, plot a histogram of the `atomization_energy` array
*    In the web-interface, answer the following question:  *What is the value range of the atomization energies?*

**Note:** The atomization energies are in units of kcal/mol.

```
import matplotlib.pyplot as plt

# Plot the histogram of the atomization energies

```

# Exercises 2.2: Training/Validation/Test split:

Before we do any training and optimization, it is a good idea to split the data in to a training set, a validation set, and a test set.

Since our representations (the Coulomb Matrix featuers, and the Bag-of-Bonds features) and the atomization energy labels are in Numpy format, you can do the splitting utilizing Numpy's "slice" notation.

Below are some examples of how to get certain rows out of a matrix which you can use for inspiration

```
# Get rows 0-9:
rows1 = my_matrix[:10]

# Get rows 10-19
rows2 = my_matrix[10:20]

# Get rows 20 and onwards
rows2 = my_matrix[20:]
```

### Question 2.2.1:  (6.2)

In this question, you have to split the three Numpy arrays `cmat_features`, `bob_features`, and `atomization_energy` each into a training, a validation and a test part. In total you will end up with 9 numpy arrays in total.

**Note:** It is important that no molecule is in more than one set.

**To-do:**
*    Split the data into Training/Validation/Test sets of roghly 60%/20%/20% of the total dataset each.

*    How many molecules are there the training, validation and test sets with a 60%/20%/20% split? In the web-interface select the right answer.
 * Since the sizes of splits can vary slightly depending on how you did the splitting, choose the option that is closest to what you found.

```
# Define the splits here

# Split the coulomb matrix features
cmat_features_training = ???
cmat_features_validation = ???
cmat_features_test = ???

# Split the bag-of-bond features
bob_features_training = ???
bob_features_validation = ???
bob_features_test = ???

# Split the atomization energies
atomization_energy_training = ???
atomization_energy_validation = ???
atomization_energy_test = ???

# Print sizes of the three data splits)
print(len(atomization_energy_training))
print(len(atomization_energy_validation))
print(len(atomization_energy_test))
```

## Exercise 2.3: Kernel Ridge Regression Model

Now that we have defined our training, validation and test splits, we are now ready to fit a kernel ridge regresssion model.

Just like the previous examples with classifiers and linear regression, Scikit-Learn has built-in support for kernel ridge regresssion. Below is some boiler-plate code you can use:

```
# Import the Machine
from sklearn.kernel_ridge import KernelRidge

# Make a machine, see text for explanation of key-words
machine = KernelRidge(alpha=1e-9, kernel="rbf", gamma=1e-4)

# Fit the machine using the training features and training labels
machine.fit(features_training, labels_training)
```
With the fitted macine, you can now make predictions for a set of test features and get the predicted labels from your machine:

```
# Predict y-values using features_training features
labels_test_predicted = machine.predict(features_test)
```

If you paid attention to the keyword arguments for the machine `KernelRidge(alpha=1e-9, kernel="rbf", gamma=1e-4)`, an explanation is given here:

1.   `alpha=1e-9` is the regularizer, i.e. the small number to add to the diagonal of the kernel matrix. $10^{-9}$ is a good value for many problems and we don't have to optimize this.
2.   `kernel="rbf"` Another word for Gaussian kernel function is also the radial basis function (RBF). This keyword tells the machine to use the Gaussian kernel.
3.   `gamma=1e-4` Gamma (i.e. $\gamma$) is the length scale of Gaussian/radial basis function.

The kernel used in Scikit-learn is defined as follows:

\begin{equation}
K\left(\mathbf{x}_i, \mathbf{x}_j\right) = \exp\left( -\gamma \|x - y \|^2 \right)
\end{equation}
Later, we will optimize the value of $\gamma$ to ensure the most accurate machine learning mode.

The accuracy of a machine learning model can be assessed, for example, by calculating the mean-absolute-error (MAE) for a test or validation set:
\begin{equation}
\text{MAE} = \frac{1}{N}\sum_{i=1}^N |y_i^{true} - y_i^{predicted} |
\end{equation}

### Question 2.3.1:  (6.3)
Make a kernel ridge regression machine as described above.

Use the same parameters values as in the example above, that is `alpha=1e-9, kernel="rbf", gamma=1e-4`.

**To-do:**
* First, train your model on the training set with Coulomb Matrix features (i.e. `cmat_features_training` and `atomization_energy_training`).
* Next predict the atomization energy on the test set (`cmat_features_test`), and calculate the mean-absolute-error (MAE) between the predicted atomization energies for the test set and the true atomization energies (stored in `atomization_energy_test`).
* In the web interface, enter the MAE of the atomization energies you calculated using Coulomb matrix features.
 * You are allowed a margin of +/- 4 kcal/mol since the numbers might vary slightly.

```
from sklearn.kernel_ridge import KernelRidge

# Implement kernel ridge regression model
# Train the model on cmat_features_training and atomization_energy_training

# Predict atomization energies on cmat_features_test

# Calculate the MAE between the predicted and true atomization energies for the test set

```

### Question 2.3.2:   (6.4)

Repeat what you did in Question 2.3.1, but this time train on the training set with Bag-of-Bonds features

* Train your kernel ridge regression using `bob_features_training` and `atomization_energy_training`.
* Next predict the atomization energy on the test set (`bob_features_test`), and calculate the mean-absolute-error (MAE) between the predicted atomization energies for the test set and the true atomization energies (stored in `atomization_energy_test`).
* In the web interface, enter the MAE of the atomization energies you calculated using Bag-of-Bonds features.
 * You are allowed a margin of +/- 4 kcal/mol since the numbers might vary slightly.

```
from sklearn.kernel_ridge import KernelRidge

# Implement kernel ridge regression model
# Train the model on bob_features_training and atomization_energy_training

# Predict atomization energies on bob_features_test

# Calculate the MAE between the predicted and true atomization energies for the test set

```

## Exercise 2.4: Hyperparameter Optimizaton.

In order to ensure the most accurate predictions we need to optimize the hyperparameters of the model. Kernel ridge regression models can be especially sensitive to the length-scale parameter, in our case the parameter $\gamma$ (`gamma`).

To set the value of gamma in the machine learning mode to, for example, $10^{-8}$ you can change the argument `gamma=` as follows:

```
machine = KernelRidge(alpha=1e-9, kernel="rbf", gamma=1e-8)
```

One common strategy to optimize such parameter is to scan all parameters over a logarithmically-space range.

To avoid overfitting (i.e. fitting parameters to the test set), it is common practice to train the model on the test set and predict on the *validation* set. The best/optimal value is the parameter with the lowest prediction error on the validation set.

### Question 2.4.1: (6.5)
Previously, we only trained machines for only for `gamma=1e-4`. In this question, you will train machines for $\gamma$ in a range between $1.0$ down to $10^{-9}$ (see code). First, we do this for the Coulomb matrix features.

**To-do:**
* For each value of `gamma`, train a machines on the training set using the Coulomb Matrix features, i.e. `cmat_features_trainig`.
* For each value of `gamma`, predict the atomization energy for the validation set, using Coulomb Matrix features, i.e. `cmat_features_validation`.
* Next, calculate the mean-absolute-error between the predicted and true atomization energies for the validation set.
* In the web-interface, select which values of `gamma` you found to give the lowest MAE using Coulomb matrix features.
 * Note: As some values of `gamma` will give MAE values that are very close, there are several right answers to this question

```
# The list of Gamma that we wish to try
gammas = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9]

# Kernel ridge regression models for each gamma on cmat_features_trainig
# Predict on the atomization energy for cmat_features_validation
# Calculate the MAE between the true and predicted atomization energies

```

### Question 2.4.2:  (6.6)
In the previous question, you found the best value of `gamma` for Coulomb matrix features. In this question you will repeat the same process and find the opimal value of `gamma` for the Bag-of-Bonds features.

**To-do:**
* For each value of `gamma`, train a machines on the training set using the Bag-of-Bonds features, i.e. `bob_features_trainig`.
* For each value of `gamma`, predict the atomization energy for the validation set, using Bag-of-Bonds features, i.e. `bob_features_validation`.
* Next, calculate the mean-absolute-error between the predicted and true atomization energies for the validation set.
* In the web-interface, select which values of `gamma` you found to give the lowest MAE using Bag-of-Bonds features.
 * Note: As some values of `gamma` will give MAE values that are very close, there are several right answers to this question

```
# The list of Gamma that we wish to try
gammas = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9]

# Kernel ridge regression models for each gamma on bob_features_trainig
# Predict on the atomization energy for bob_features_validation
# Calculate the MAE between the true and predicted atomization energies

```

## Exercise 2.5: Learning Curves

In order to see which of the two reprsentations  (Coulomb Matrix and Bag-of-Bonds) gives the most accurate kernel ridge regression model, we can compare learning curves for the machines based on the two representations.

### Question 2.5.1: (6.7)
For both representations, train machines using 500, 1000, 2000, and 4000 molecules from the training set. Use the best value of `gamma` which you found in the previous exercise for each machine.

Again, use Numpy's slice notation to extract only the neccesary number of rows from the Numpy arrays in the training set.

**To-do:**
*   Train machines for each training set size with the two representation (Coulomb Matrix and Bag-of-bonds)
*   For each machine you train, predict the atomization energy for the test set and calculate the mean-absolue-error (MAE) to the true energies for the test set.
*   Plot the MAE values as a function of the training set size.
*   The representation which yields the most accurate predictions is the one with the lowest learning curve. In the web interface, select the representation that gave you the learning curves with the lowest MAE values.

```
# Implement learning curves for kernel ridge regression here

```

ext: ipynb

max_stars_repo_name: andersx/python-intro

max_issues_repo_name: andersx/ml-intro

max_forks_repo_name: andersx/ml-intro


text:
#  CS109B Data Science 2: Advanced Topics in Data Science

## Lab 4 - Bayesian Analysis

**Harvard University**

**Spring 2020**

**Instructors:** Mark Glickman, Pavlos Protopapas, and Chris Tanner

**Lab Instructors:** Chris Tanner and Eleni Angelaki Kaxiras

**Content:** Eleni Angelaki Kaxiras

---

```python
## RUN THIS CELL TO PROPERLY HIGHLIGHT THE EXERCISES
import requests
from IPython.core.display import HTML
styles = requests.get("https://raw.githubusercontent.com/Harvard-IACS/2019-CS109B/master/content/styles/cs109.css").text
HTML(styles)
```

```python
import pymc3 as pm
from pymc3 import summary
```

```python
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
import pandas as pd
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')
```

```python
print('Running on PyMC3 v{}'.format(pm.__version__))
```

    Running on PyMC3 v3.8

```javascript
%%javascript
IPython.OutputArea.auto_scroll_threshold = 20000;
```

    <IPython.core.display.Javascript object>

## Learning Objectives

By the end of this lab, you should be able to:
* Understand how probability distributions work.
* Apply Bayes Rule in calculating probabilities.
* Understand how to apply Bayesian analysis using PyMC3
* Avoid getting fired when talking to your Bayesian employer.

**This lab corresponds to Lectures 6, 7, and 8, and maps to Homework 3.**

## Table of Contents

1. The Bayesian Way of Thinking or Is this a Fair Coin?
2. [Intro to `pyMC3`](#pymc3).
3. [Bayesian Linear Regression](#blr).
4. [Try this at Home: Example on Mining Disasters](#no4).

## 1. The Bayesian way of Thinking

```
Here is my state of knowledge about the situation. Here is some data, I am now going to revise my state of knowledge.
```

Table Exercise: Discuss the statement above with your table mates and make sure everyone understands what it means and what constitutes Bayesian way of thinking. Finally, count the Bayesians among you.

### A. Bayes Rule

\begin{equation}
\label{eq:bayes}
P(A|\textbf{B}) = \frac{P(\textbf{B} |A) P(A) }{P(\textbf{B})}
\end{equation}

$P(A|\textbf{B})$ is the **posterior** distribution, prob(hypothesis | data)

$P(\textbf{B} |A)$ is the **likelihood** function, how probable is my data **B** for different values of the parameters

$P(A)$ is the marginal probability to observe the data, called the **prior**, this captures our belief about the data before observing it.

$P(\textbf{B})$ is the marginal distribution (sometimes called marginal likelihood)

Table Exercise: Solve the Monty Hall Paradox using Bayes Rule.

You are invited to play a game. There are 3 doors behind **one** of which are the keys to a brand new red Tesla. There is a goat behind each of the other two.

You are asked to pick one door, and let's say you pick **Door1**. The host knows where the keys are. Of the two remaining closed doors, he will always open the door that has a goat behind it. He'll say "I will do you a favor and open **Door2**". So he opens Door2 inside which there is, of course, a goat. He now asks you, do you want to open the initial Door you chose or change to **Door3**? Generally, in this game, when you are presented with this choice should you swap the doors?

**Initial Steps:**
- Start by defining the `events` of this probabilities game. One definition is:

    - $A_i$: car is behind door $i$

    - $B_i$ host opens door $i$

$i\in[1,2,3]$

- In more math terms, the question is: is the probability that the price is behind **Door 1** higher than the probability that the price is behind **Door2**, given that an event **has occured**?

### B. Bayes Rule written with Probability Distributions

We have data that we believe come from an underlying distribution of unknown parameters. If we find those parameters, we know everything about the process that generated this data and we can make inferences (create new data).

$$
\begin{equation}
\label{eq:bayes}
P(\theta|\textbf{D}) = \frac{P(\textbf{D} |\theta) P(\theta) }{P(\textbf{D})}
\end{equation}
$$

#### But what is $\theta \;$?

$\theta$ is an unknown yet fixed set of parameters. In Bayesian inference we express our belief about what $\theta$ might be and instead of trying to guess $\theta$ exactly, we look for its **probability distribution**. What that means is that we are looking for the **parameters** of that distribution. For example, for a Poisson distribution our $\theta$ is only $\lambda$. In a normal distribution, our $\theta$ is often just $\mu$ and $\sigma$.

### C. A review of Common Probability Distributions

#### Discrete Distributions

The random variable has a **probability mass function (pmf)** which measures the probability that our random variable will take a specific value $y$, denoted $P(Y=y)$.

- **Bernoulli** (binary outcome, success has probability $\theta$, $one$ trial):
$
P(Y=k) =  \theta^k(1-\theta)^{1-k}
$

- **Binomial** (binary outcome, success has probability $\theta$, $n$ trials):
\begin{equation}
P(Y=k) =  {{n}\choose{k}} \cdot \theta^k(1-\theta)^{n-k}
\end{equation}

*Note*: Binomial(1,$p$) = Bernouli($p$)

- **Negative Binomial**

- **Poisson** (counts independent events occurring at a rate)
\begin{equation}
P\left( Y=y|\lambda \right) = \frac{{e^{ - \lambda } \lambda ^y }}{{y!}}
\end{equation}
y = 0,1,2,...

- **Discrete Uniform**

- **Categorical, or Multinulli** (random variables can take any of K possible categories, each having its own probability; this is a generalization of the Bernoulli distribution for a discrete variable with more than two possible outcomes, such as the roll of a die)

- **Dirichlet-multinomial** (a generalization of the beta distribution for many variables)

#### Continuous Distributions

The random variable has a **probability density function (pdf)**.
- **Uniform** (variable equally likely to be near each value in interval $(a,b)$)
\begin{equation}
P(X = x) = \frac{1}{b - a}
\end{equation}
anywhere within the interval $(a, b)$, and zero elsewhere.

- **Normal** (a.k.a. Gaussian)
\begin{equation}
X \sim  \mathcal{N}(\mu,\,\sigma^{2})
\end{equation}

    A Normal distribution can be parameterized either in terms of precision $\tau$ or standard deviation ($\sigma^{2}$. The link between the two is given by
\begin{equation}
\tau = \frac{1}{\sigma^{2}}
\end{equation}
 - Mean $\mu$
 - Variance $\frac{1}{\tau}$ or $\sigma^{2}$
 - Parameters: `mu: float`, `sigma: float` or `tau: float`

- **Beta** (variable ($\theta$) taking on values in the interval $[0,1]$, and parametrized by two positive parameters, $\alpha$ and $\beta$ that control the shape of the distribution.

*Note:*Beta is a good distribution to use for priors (beliefs) because its range is $[0,1]$ which is the natural range for a probability and because we can model a wide range of functions by changing the $\alpha$ and $\beta$ parameters.

\begin{equation}
\label{eq:beta}
P(\theta) = \frac{1}{B(\alpha, \beta)} {\theta}^{\alpha - 1} (1 - \theta)^{\beta - 1} \propto {\theta}^{\alpha - 1} (1 - \theta)^{\beta - 1}
\end{equation}

where the normalisation constant, $B$, is a beta function of $\alpha$ and $\beta$,

\begin{equation}
B(\alpha, \beta) = \int_{t=0}^1 t^{\alpha - 1} (1 - t)^{\beta - 1} dt.
\end{equation}

- **Exponential**

- **Gamma**

 #### Code Resources:
 - Statistical Distributions in numpy/scipy: [scipy.stats](https://docs.scipy.org/doc/scipy/reference/stats.html)
 - Statistical Distributions in pyMC3: [distributions in PyMC3](https://docs.pymc.io/api/distributions.html) (we will see those below).

Exercise: Plot a Discrete variable

Change the value of $\mu$ in the Poisson PMF and see how the plot changes. Remember that the y-axis in a discrete probability distribution shows the probability of the random variable having a specific value in the x-axis.

\begin{equation}
P\left( X=k \right) = \frac{{e^{ - \mu } \mu ^k }}{{k!}}
\end{equation}

**stats.poisson.pmf(x, mu)** $\mu$(mu) is our $\theta$ in this case.

```python
plt.style.use('seaborn-darkgrid')
x = np.arange(0, 30)
for m in [0.5, 3, 8]:
    pmf = stats.poisson.pmf(x, m)
    plt.plot(x, pmf, 'o', alpha=0.5, label='$\mu$ = {}'.format(m))
plt.xlabel('random variable', fontsize=12)
plt.ylabel('probability', fontsize=12)
plt.legend(loc=1)
plt.ylim=(-0.1)
plt.show()
```

```python
# same for binomial
plt.style.use('seaborn-darkgrid')
x = np.arange(0, 22)
ns = [10, 17]
ps = [0.5, 0.7]
for n, p in zip(ns, ps):
    pmf = stats.binom.pmf(x, n, p)
    plt.plot(x, pmf, 'o', alpha=0.5, label='n = {}, p = {}'.format(n, p))
plt.xlabel('x', fontsize=14)
plt.ylabel('f(x)', fontsize=14)
plt.legend(loc=1)
plt.show()
```

```python
# discrete uniform
plt.style.use('seaborn-darkgrid')
ls = [0]
us = [3] # watch out, this number can only be integer!
for l, u in zip(ls, us):
    x = np.arange(l, u+1)
    pmf = [1.0 / (u - l + 1)] * len(x)
    plt.plot(x, pmf, '-o', label='lower = {}, upper = {}'.format(l, u))
plt.xlabel('x', fontsize=12)
plt.ylabel('probability P(x)', fontsize=12)
plt.legend(loc=1)
plt.show()
```

Exercise: Plot a continuous variable

Change the value of $\mu$ in the Uniform PDF and see how the plot changes.

Remember that the y-axis in a continuous probability distribution does not shows the actual probability of the random variable having a specific value in the x-axis because that probability is zero!. Instead, to see the probability that the variable is within a small margin we look at the integral below the curve of the PDF.

The uniform is often used as a noninformative prior.

```
Uniform - numpy.random.uniform(a=0.0, b=1.0, size)
```

$\alpha$ and $\beta$ are our parameters. `size` is how many tries to perform.
Our $\theta$ is basically the combination of the parameters a,b. We can also call it
\begin{equation}
\mu = (a+b)/2
\end{equation}

```python
uniform.pdf([.1,2])
```

```python
from scipy.stats import uniform

r = uniform.rvs(size=1000)
plt.plot(r, uniform.pdf(r),'r-', lw=5, alpha=0.6, label='uniform pdf')
plt.hist(r, density=True, histtype='stepfilled', alpha=0.2)
plt.ylabel(r'probability density')
plt.xlabel(f'random variable')
plt.legend(loc='best', frameon=False)
plt.show()
```

```python
from scipy.stats import beta

alphas = [0.5, 1.5, 3.0]
betas = [0.5, 1.5, 3.0]
x = np.linspace(0, 1, 1000)
colors = ['red', 'green', 'blue']

fig, ax = plt.subplots(figsize=(8, 5))

for a, b, colors in zip(alphas, betas, colors):
    dist = beta(a, b)
    plt.plot(x, dist.pdf(x), c=colors,
             label=f'a={a}, b={b}')

ax.set_ylim(0, 3)

ax.set_xlabel(r'$\theta$')
ax.set_ylabel(r'$p(\theta|\alpha,\beta)$')
ax.set_title('Beta Distribution')

ax.legend(loc='best')
fig.show();
```

```python
plt.style.use('seaborn-darkgrid')
x = np.linspace(-5, 5, 1000)
mus = [0., 0., 0., -2.]
sigmas = [0.4, 1., 2., 0.4]
for mu, sigma in zip(mus, sigmas):
    pdf = stats.norm.pdf(x, mu, sigma)
    plt.plot(x, pdf, label=r'$\mu$ = '+ f'{mu},' + r'$\sigma$ = ' + f'{sigma}')
plt.xlabel('random variable', fontsize=12)
plt.ylabel('probability density', fontsize=12)
plt.legend(loc=1)
plt.show()
```

```python
plt.style.use('seaborn-darkgrid')
x = np.linspace(-5, 5, 1000)
mus = [0., 0., 0., -2.] # mean
sigmas = [0.4, 1., 2., 0.4] # std
for mu, sigma in zip(mus, sigmas):
      plt.plot(x, uniform.pdf(x, mu, sigma), lw=5, alpha=0.4, \
               label=r'$\mu$ = '+ f'{mu},' + r'$\sigma$ = ' + f'{sigma}')
plt.xlabel('random variable', fontsize=12)
plt.ylabel('probability density', fontsize=12)
plt.legend(loc=1)
plt.show()
```

### D. Is this a Fair Coin?

We do not want to promote gambling but let's say you visit the casino in **Monte Carlo**. You want to test your theory that casinos are dubious places where coins have been manipulated to have a larger probability for tails. So you will try to estimate how fair a coin is based on 100 flips.

You begin by flipping the coin. You get either Heads ($H$) or Tails ($T$) as our observed data and want to see if your posterior probabilities change as you obtain more data, that is, more coin flips. A nice way to visualize this is to plot the posterior probabilities as we observe more flips (data).

We will be using Bayes rule. $\textbf{D}$ is our data.

$
\begin{equation}
\label{eq:bayes}
P(\theta|\textbf{D}) = \frac{P(\textbf{D} |\theta) P(\theta) }{P(\textbf{D})}
\end{equation}
$

In the case of a coin toss when we observe $k$ heads in $n$ tosses:

$
\begin{equation}
\label{eq:bayes}
P(\theta|\textbf{k}) = Beta(\alpha + \textbf{k}, \beta + n - \textbf{k})
\end{equation}
$

we can say that $\alpha$ and $\beta$ play the roles of a "prior number of heads" and "prior number of tails".

```python
# play with the priors - here we manually set them but we could be sampling from a separate Beta
trials = np.array([0, 1, 3, 5, 10, 15, 20, 100, 200, 300])
heads = np.array([0, 1, 2, 4, 8, 10, 10, 50, 180, 150])
x = np.linspace(0, 1, 100)

# for simplicity we set a,b=1

plt.figure(figsize=(10,8))
for k, N in enumerate(trials):
    sx = plt.subplot(len(trials)/2, 2, k+1)
    posterior = stats.beta.pdf(x, 1 + heads[k], 1 + trials[k] - heads[k])
    plt.plot(x, posterior, alpha = 0.5, label=f'{trials[k]} tosses\n {heads[k]} heads');
    plt.fill_between(x, 0, posterior, color="#348ABD", alpha=0.4)
    plt.legend(loc='upper left', fontsize=10)
    plt.legend()
    plt.autoscale(tight=True)

plt.suptitle("Posterior probabilities for coin flips", fontsize=15);
plt.tight_layout()
plt.subplots_adjust(top=0.88)
```

 [Top](#top)

## 2. Introduction to `pyMC3`

PyMC3 is a Python library for programming Bayesian analysis, and more specifically, data creation, model definition, model fitting, and posterior analysis. It uses the concept of a `model` which contains assigned parametric statistical distributions to unknown quantities in the model. Within models we define random variables and their distributions. A distribution requires at least a `name` argument, and other `parameters` that define it. You may also use the `logp()` method in the model to build the model log-likelihood function. We define and fit the model.

PyMC3 includes a comprehensive set of pre-defined statistical distributions that can be used as model building blocks. Although they are not meant to be used outside of a `model`, you can invoke them by using the prefix `pm`, as in `pm.Normal`.

#### Markov Chain Monte Carlo (MCMC) Simulations

PyMC3 uses the **No-U-Turn Sampler (NUTS)** and the **Random Walk Metropolis**, two Markov chain Monte Carlo (MCMC) algorithms for sampling in posterior space. Monte Carlo gets into the name because when we sample in posterior space, we choose our next move via a pseudo-random process. NUTS is a sophisticated algorithm that can handle a large number of unknown (albeit continuous) variables.

```python
with pm.Model() as model:
    z = pm.Normal('z', mu=0., sigma=5.)
    x = pm.Normal('x', mu=z, sigma=1., observed=5.)
print(x.logp({'z': 2.5}))
print(z.random(10, 100)[:10])
```

**References**:

- *Salvatier J, Wiecki TV, Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:e55* [(https://doi.org/10.7717/peerj-cs.55)](https://doi.org/10.7717/peerj-cs.55)
- [Distributions in PyMC3](https://docs.pymc.io/api/distributions.html)
- [More Details on Distributions](https://docs.pymc.io/developer_guide.html)

Information about PyMC3 functions including descriptions of distributions, sampling methods, and other functions, is available via the `help` command.

```python
#help(pm.Poisson)
```

 [Top](#top)

## 3. Bayesian Linear Regression

Let's say we want to predict outcomes Y as normally distributed observations with an expected value $mu$ that is a linear function of two predictor variables, $\bf{x}_1$ and $\bf{x}_2$.

\begin{equation}
\mu = \alpha + \beta_1 \bf{x}_1 + \beta_2 x_2
\end{equation}

\begin{equation}
Y \sim  \mathcal{N}(\mu,\,\sigma^{2})
\end{equation}

where $\sigma^2$ represents the measurement error.

In this example, we will use $\sigma^2 = 10$

We also choose the parameters as normal distributions:

\begin{eqnarray}
\alpha \sim  \mathcal{N}(0,\,10) \\
\beta_i \sim  \mathcal{N}(0,\,10) \\
\sigma^2 \sim  |\mathcal{N}(0,\,10)|
\end{eqnarray}

We will artificially create the data to predict on. We will then see if our model predicts them correctly.

```python
# Initialize random number generator
np.random.seed(123)

# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]

# Size of dataset
size = 100

# Predictor variable
X1 = np.linspace(0, 1, size)
X2 = np.linspace(0,.2, size)

# Simulate outcome variable
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma

fig, ax = plt.subplots(1,2, figsize=(10,6), sharex=True)
ax[0].scatter(X1,Y)
ax[1].scatter(X2,Y)
ax[0].set_xlabel(r'$x_1$', fontsize=14)
ax[0].set_ylabel(r'$Y$', fontsize=14)
ax[1].set_xlabel(r'$x_2$', fontsize=14)
ax[1].set_ylabel(r'$Y$', fontsize=14)
```

```python
from pymc3 import Model, Normal, HalfNormal

basic_model = Model()

with basic_model:

    # Priors for unknown model parameters, specifically create stochastic random variables
    # with Normal prior distributions for the regression coefficients,
    # and a half-normal distribution for the standard deviation of the observations, σ.
    alpha = Normal('alpha', mu=0, sd=10)
    beta = Normal('beta', mu=0, sd=10, shape=2)
    sigma = HalfNormal('sigma', sd=1)

    # Expected value of outcome - posterior
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood (sampling distribution) of observations
    Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)
```

```python
# model fitting with sampling
from pymc3 import NUTS, sample, find_MAP
from scipy import optimize

with basic_model:

    # obtain starting values via MAP
    start = find_MAP(fmin=optimize.fmin_powell)

    # instantiate sampler
    step = NUTS(scaling=start)

    # draw 2000 posterior samples
    trace = sample(2000, step, start=start)
```

```python
from pymc3 import traceplot

traceplot(trace);
```

```python
results = pm.summary(trace,
                     var_names=['alpha', 'beta', 'sigma'])
results
```

This linear regression example is from the original paper on PyMC3: *Salvatier J, Wiecki TV, Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2:e55 https://doi.org/10.7717/peerj-cs.55*

 [Top](#top)

## 4. Try this at Home: Example on Mining Disasters
We will go over the classical `mining disasters from 1851 to 1962` dataset.

This example is from the [pyMC3 Docs](https://docs.pymc.io/notebooks/getting_started.html).

```python
import pandas as pd
disaster_data = pd.Series([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
                           3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
                           2, 2, 3, 4, 2, 1, 3, np.nan, 2, 1, 1, 1, 1, 3, 0, 0,
                           1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
                           0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
                           3, 3, 1, np.nan, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
                           0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
fontsize = 12
years = np.arange(1851, 1962)
plt.figure(figsize=(10,5))
#plt.scatter(years, disaster_data);
plt.bar(years, disaster_data)
plt.ylabel('Disaster count', size=fontsize)
plt.xlabel('Year', size=fontsize);
plt.title('Was there a Turning Point in Mining disasters from 1851 to 1962?', size=15);
```

#### Building the model

**Step1:** We choose the probability model for our experiment. Occurrences of disasters in the time series is thought to follow a **Poisson** process with a large **rate** parameter in the early part of the time series, and from one with a smaller **rate** in the later part. We are interested in locating the change point in the series, which perhaps is related to changes in mining safety regulations.

```
disasters = pm.Poisson('disasters', rate, observed=disaster_data)
```

We have two rates, `early_rate` if $t<=s$, and `late_rate` if $t>s$, where $s$ is the year the switch was made (a.k.a. the `switchpoint`).

**Step2:** Choose a prior distributions of the two rates, what we believe the rates were before we observed the data, and the switchpoint. We choose Exponential.
```
early_rate = pm.Exponential('early_rate', 1)
```

The parameters of this model are:

**Note:** Watch for missing values. Missing values are handled transparently by passing a MaskedArray or a pandas.DataFrame. Behind the scenes, another random variable, disasters.missing_values is created to model the missing values. If you pass a np.array with missing values you will get an error.

```python
with pm.Model() as disaster_model:

    # discrete
    switchpoint = pm.DiscreteUniform('switchpoint', lower=years.min(), upper=years.max(), testval=1900)

    # Priors for pre- and post-switch rates number of disasters
    early_rate = pm.Exponential('early_rate', 1)
    late_rate = pm.Exponential('late_rate', 1)

    # our theta - allocate appropriate Poisson rates to years before and after current
    # switch is an `if` statement in puMC3
    rate = pm.math.switch(switchpoint >= years, early_rate, late_rate)

    # our observed data as a likelihood function of the `rate` parameters
    # shows how we think our data is distributed
    disasters = pm.Poisson('disasters', rate, observed=disaster_data)
```

#### Model Fitting

```python
# there are defaults but we can also more explicitly set the sampling algorithms
with disaster_model:

    # for continuous variables
    step1 = pm.NUTS([early_rate, late_rate])

    # for discrete variables
    step2 = pm.Metropolis([switchpoint, disasters.missing_values[0]] )

    trace = pm.sample(10000, step=[step1, step2])
    # try different number of samples
    #trace = pm.sample(5000, step=[step1, step2])
```

    Multiprocess sampling (4 chains in 4 jobs)
    CompoundStep
    >NUTS: [late_rate, early_rate]
    >CompoundStep
    >>Metropolis: [disasters_missing]
    >>Metropolis: [switchpoint]
    Sampling 4 chains, 0 divergences: 100%|██████████| 42000/42000 [00:13<00:00, 3189.17draws/s]
    The number of effective samples is smaller than 10% for some parameters.

#### Posterior Analysis

On the left side plots we notice that our early rate is between 2.5 and 3.5 disasters a year. In the late period it seems to be between 0.6 and 1.2 so definitely lower.

The right side plots show the samples we drew to come to our conclusion.

```python
pm.traceplot(trace, ['early_rate', 'late_rate', 'switchpoint'], figsize=(20,10));
```

```python
results = pm.summary(trace,
                     var_names=['early_rate', 'late_rate', 'switchpoint'])
results
```

      mean
      sd
      hpd_3%
      hpd_97%
      mcse_mean
      mcse_sd
      ess_mean
      ess_sd
      ess_bulk
      ess_tail
      r_hat

      early_rate
      3.089
      0.290
      2.540
      3.624
      0.002
      0.001
      19436.0
      19436.0
      19343.0
      23587.0
      1.0

      late_rate
      0.932
      0.119
      0.715
      1.156
      0.001
      0.001
      19680.0
      19680.0
      19355.0
      24082.0
      1.0

      switchpoint
      1889.730
      2.426
      1886.000
      1894.000
      0.057
      0.040
      1807.0
      1807.0
      1832.0
      3126.0
      1.0

```python

```


text:
# Bouncing ball with the Euler integrator
As an introduction to solving ODEs we are solving the problem of the bouncing ball with the simplest of all integrators: the forward Euler scheme.

We have to solve the *second order ODE* (Newton's equations of motion with constant acceleration)
$$
\frac{d^2 y}{dt^2} = -g.
$$

The Euler scheme for any *first order ODE*

$$
\frac{dy}{dt} = f(y, t)
$$

is
$$
y(t + h) = y(t) + h f(y(t), t).
$$

In order to solve the original 2nd order equation of motion we make use of the fact that one $n$-th order ODE can be written as $n$ coupled first order ODEs, namely

\begin{align}
\frac{dy}{dt} &= v\\
\frac{dv}{dt} &= -g.
\end{align}

Solve each of the first order ODEs with Euler:

\begin{align}
y(t + h) &= y(t) + h v(t)\\
v(t + h) &= v(t) - h g.
\end{align}

```python
import numpy as np
```

## Free fall

Start with free fall

```python
g = -9.81
y = 0.0
v = 0.0

t = 0
dt = 0.01

data = []

while t < 10:
    y = y + v*dt
    v = v + g*dt
    data.append([t, y, v])
    t += dt

data = np.array(data)
```

```python
data.shape
```

    (1001, 3)

```python
data = data.transpose()
data.shape  # t, y, v
```

    (3, 1001)

```python
import matplotlib.pyplot as plt
%matplotlib inline
```

```python
plt.plot(data[0], data[1])
plt.xlabel("time (s)")
plt.ylabel("position (m)")
```

## Bouncing
Add a floor at $y = -5$.

What happens at the floor? – The velocity changes (elastic collision).

```python
g = -9.81
y = 0.0
v = 0.0

t = 0
dt = 0.01

y_floor = -5

data = []

while t < 10:
    y += v*dt
    if y > y_floor:
        v += g*dt
    else:
        v = -v   # bounce off floor
    data.append([t, y, v])
    t += dt

data = np.array(data).transpose()
```

```python
plt.plot(data[0], data[1])
plt.xlabel("time (s)")
plt.ylabel("position (m)")
```

```python

```