<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Computational Bayes Final Project
Nathan Yee
Uma Desai  
To understand Bayesian Linear Regression, we read through and used many ideas from Cypress Frankenfeld. http
Step1:  Let's start by loading our data into a Pandas dataframe to see what we're working with.
Step2:  Next, let's create vectors of our ages and heights.
Step3:  Now let's visualize our data to make sure that linear regression is appropriate for predicting its distributions.
Step5:  Our data looks pretty linear. We can now calculate the slope and intercept of the line of least squares. We abstract numpy's least squares function using a function of our own.
Step6:  To use our leastSquares function, we input our age and height vectors as our x and y arguments. Next, let's call leastSquares to get the slope and intercept, and use the slope and intercept to calculate the size of our intercept and slope ranges.
Step7:  Now, we can visualize the slope and intercept on the same plot as the data to make sure it is working correctly.
Step8:  Looks great! Based on the plot above, we are confident that Bayesian linear regression will give us reasonable distributions for predicting future values. Now we need to create our hypotheses. Each hypothesis will consist of a range of intercepts, slopes and sigmas.
Step9:  Loop through ages and heights to create tuples as data for our hypotheses.
Step11:  Next make a least squares class that inherits from Suite and Joint where likelihood is calculated based on error from data. The likelihood function will depend on the data and normal distributions for each hypothesis.
Step12:  Now instantiate a LeastSquaresHypos suite with our hypos.
Step13:  And update the suite with our data. We choose to update a single item at a time to normalize after every item in data.
Step14:  We can now plot marginal distributions to visualize the probability distribution for each of our hypotheses for intercept, slope, and sigma values. Our hypotheses were carefully picked based on ranges that we found worked well, which is why all the intercepts, slopes, and sigmas that are important to this dataset are included in our hypotheses.
Step15:  Since our slopes are relatively close to our original slope, we are relatively confdient that these produce reasonable distributions even though it does get cut off.
Step18:  Next, we want to sample random data from our hypotheses. To do this, we will make two functions, getHeights and getRandomData. getRandomData calls getHeights to obtain random height values.
Step19:  Now we take 100000 random samples of pairs of months and heights. Here we want at least 100000 items so that we can get very smooth sampling.
Step20:  Now plot the data to see what the randomly sampled data looks like beneath the original data and line of least squares. We notice that we have trouble seeing the exact size of our distribution as points overlap each other. To see what the distribution actually looks like, we must make a density plot.
Step21:  Next, we want to get the intensity of the data at locations. We do that by adding the randomly sampled values to buckets. This gives us intensity values for a grid of pixels in our sample range.
Step22:  Next, we unpack the buckets into three vectors, pcolor_months, pcolor_heights, and pcolor_densities
Step29:  Since density plotting is much simpler in Mathematica, we are going to export all our data to csv files and plot them in Mathematica
Step31:  Now, we have a figure that contains superimposed density plot, original data scatter plot, and least squares line plot. While the density plot is not strictly necessary, it gives us a much better idea as to how large of a distribution we have. Note that the x axis is ages (months) and the y axis is height (cm).
<img src="ageHeightAllPlots5.png" alt="Density Plot with orignal data/fit" height="400" width="400">
So, now we can go back to our original question. What is the probability that a child between the ages of 18 - 29 is less than 76.1 cm (the lowest height in our dataset)? To answer this question, we sample our random data and see whether or not it is less than 76.1 cm. We will do this using a function.
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function, division
% matplotlib inline
import warnings
warnings.filterwarnings('ignore')
import math
import numpy as np
from thinkbayes2 import Pmf, Cdf, Suite, Joint, EvalNormalPdf
import thinkplot
import pandas as pd
import matplotlib.pyplot as plt
Explanation: Computational Bayes Final Project
Nathan Yee
Uma Desai  
To understand Bayesian Linear Regression, we read through and used many ideas from Cypress Frankenfeld. http://allendowney.blogspot.com/2015/04/two-hour-marathon-by-2041-probably.html
Height of Children in Kalama, Egypt
From: http://lib.stat.cmu.edu/DASL/Datafiles/Ageandheight.html
The height of a child is not stable but increases over time. Since the pattern of growth varies from child to child, one way to understand the general growth pattern is by using the average of several children's heights, as presented in this data set. The scatterplot of height versus age is almost a straight line, showing a linear growth pattern. The straightforward relationship between height and age provides a simple illustration of linear relationships, correlation, and simple regression. 
Description: Mean heights of a group of children in Kalama, an Egyptian village that is the site of a study of nutrition in developing countries. The data were obtained by measuring the heights of all 161 children in the village each month over several years. 
Age: Age in months
Height: Mean height in centimeters for children at this age 
What is the probability that a child between the ages of 18 - 29 is less than 76.1 cm (the lowest height in our dataset)?
End of explanation
df = pd.read_csv('ageVsHeight.csv', skiprows=0, delimiter='\t')
df
Explanation: Let's start by loading our data into a Pandas dataframe to see what we're working with.
End of explanation
ages = np.array(df['age'])
heights = np.array(df['height'])
Explanation: Next, let's create vectors of our ages and heights.
End of explanation
plt.plot(ages, heights, 'o', label='Original data', markersize=10)
plt.xlabel('Age (months)')
plt.ylabel('Height (cm)')
Explanation: Now let's visualize our data to make sure that linear regression is appropriate for predicting its distributions.
End of explanation
def leastSquares(x, y):
    
    leastSquares takes in two arrays of values. Then it returns the slope and intercept
    of the least squares of the two.
    
    Args:
        x (numpy array): numpy array of values.
        y (numpy array): numpy array of values.
        
    Returns:
        slope, intercept (tuple): returns a tuple of floats.
    
    A = np.vstack([x, np.ones(len(x))]).T
    slope, intercept = np.linalg.lstsq(A, y)[0]
    return slope, intercept
Explanation: Our data looks pretty linear. We can now calculate the slope and intercept of the line of least squares. We abstract numpy's least squares function using a function of our own.
End of explanation
slope, intercept = leastSquares(ages, heights)
print(slope, intercept)
intercept_range = .03 * intercept
slope_range = .05 * slope
Explanation: To use our leastSquares function, we input our age and height vectors as our x and y arguments. Next, let's call leastSquares to get the slope and intercept, and use the slope and intercept to calculate the size of our intercept and slope ranges.
End of explanation
plt.plot(ages, heights, 'o', label='Original data', markersize=10)
plt.plot(ages, slope*ages + intercept, 'r', label='Fitted line')
plt.legend()
plt.xlabel('Age (months)')
plt.ylabel('Height (cm)')
plt.show()
Explanation: Now, we can visualize the slope and intercept on the same plot as the data to make sure it is working correctly.
End of explanation
intercepts = np.linspace(intercept - intercept_range, intercept + intercept_range, 20)
slopes = np.linspace(slope - slope_range, slope + slope_range, 20)
sigmas = np.linspace(2, 4, 15)
hypos = ((intercept, slope, sigma) for intercept in intercepts 
         for slope in slopes for sigma in sigmas)
Explanation: Looks great! Based on the plot above, we are confident that Bayesian linear regression will give us reasonable distributions for predicting future values. Now we need to create our hypotheses. Each hypothesis will consist of a range of intercepts, slopes and sigmas.
End of explanation
data = [(age, height) for age in ages for height in heights]
Explanation: Loop through ages and heights to create tuples as data for our hypotheses.
End of explanation
class leastSquaresHypos(Suite, Joint):
    def Likelihood(self, data, hypo):
        
        Likelihood calculates the probability of a particular line (hypo)
        based on data (ages Vs height) of our original dataset. This is
        done with a normal pmf as each hypo also contains a sigma.
        
        Args:
            data (tuple): tuple that contains ages (float), heights (float)
            hypo (tuple): intercept (float), slope (float), sigma (float)
            
        Returns:
            P(data|hypo)
        
        intercept, slope, sigma = hypo
        total_likelihood = 1
        for age, measured_height in data:
            hypothesized_height = slope * age + intercept
            error = measured_height - hypothesized_height
            total_likelihood *= EvalNormalPdf(error, mu=0, sigma=sigma)
        return total_likelihood
        
Explanation: Next make a least squares class that inherits from Suite and Joint where likelihood is calculated based on error from data. The likelihood function will depend on the data and normal distributions for each hypothesis.
End of explanation
LeastSquaresHypos = leastSquaresHypos(hypos)
Explanation: Now instantiate a LeastSquaresHypos suite with our hypos.
End of explanation
for item in data:
    LeastSquaresHypos.Update([item])
Explanation: And update the suite with our data. We choose to update a single item at a time to normalize after every item in data.
End of explanation
marginal_intercepts = LeastSquaresHypos.Marginal(0)
thinkplot.hist(marginal_intercepts)
thinkplot.Config(xlabel='intercept (cm)', ylabel='prob')
Explanation: We can now plot marginal distributions to visualize the probability distribution for each of our hypotheses for intercept, slope, and sigma values. Our hypotheses were carefully picked based on ranges that we found worked well, which is why all the intercepts, slopes, and sigmas that are important to this dataset are included in our hypotheses.
End of explanation
marginal_slopes = LeastSquaresHypos.Marginal(1)
thinkplot.hist(marginal_slopes)
thinkplot.Config(xlabel='slope (cm / month)', ylabel='prob')
marginal_sigmas = LeastSquaresHypos.Marginal(2)
thinkplot.hist(marginal_sigmas)
thinkplot.Config(xlabel='sigma', ylabel='prob')
Explanation: Since our slopes are relatively close to our original slope, we are relatively confdient that these produce reasonable distributions even though it does get cut off.
End of explanation
def getHeights(hypo_samples, random_months):
    
    getHeights takes in random hypos and random months and returns the coorisponding
    random height
    
    Args:
        hypo_samples (sequence): Sampled hypotheses from the LeastSquaredHypos suite
        random_months (sequence): Randomly sampled months
    
    random_heights = np.zeros(len(random_months))
    for i in range(len(random_heights)):
        intercept = hypo_samples[i][0]
        slope = hypo_samples[i][1]
        sigma = hypo_samples[i][2]
        month = random_months[i]
        random_heights[i] = np.random.normal((slope * month + intercept), sigma, 1)
    return random_heights
def getRandomData(start_month, end_month, n, LeastSquaresHypos):
    
    getRadomData computes and returns two sequences
    
    Args:
        start_month: (number): Starting month
        end_month: (number): Ending month
        n (int): Number of samples
        LeastSquaresHypos (Suite): Contains the hypos we want to sample
    
    Returns:
        Sequence of random_months and random_heights
    
    
    random_hypos = LeastSquaresHypos.Sample(n)
    random_months = np.random.uniform(start_month, end_month, n)
    random_heights = getHeights(random_hypos, random_months)
    return random_months, random_heights
Explanation: Next, we want to sample random data from our hypotheses. To do this, we will make two functions, getHeights and getRandomData. getRandomData calls getHeights to obtain random height values.
End of explanation
num_samples = 100000
random_months, random_heights = getRandomData(18, 29, num_samples, LeastSquaresHypos)
Explanation: Now we take 100000 random samples of pairs of months and heights. Here we want at least 100000 items so that we can get very smooth sampling.
End of explanation
plt.plot(random_months, random_heights, 'o', label='Random Sampling')
plt.plot(ages, heights, 'o', label='Original data', markersize=10)
plt.plot(ages, slope*ages + intercept, 'r', label='Fitted line')
plt.xlabel('Age (months)')
plt.ylabel('Height (cm)')
plt.legend()
plt.show()
Explanation: Now plot the data to see what the randomly sampled data looks like beneath the original data and line of least squares. We notice that we have trouble seeing the exact size of our distribution as points overlap each other. To see what the distribution actually looks like, we must make a density plot.
End of explanation
num_buckets = 70 # num_buckets^2 is actual number
# create horizontal and vertical linearly spaced ranges as buckets.
hori_range, hori_step = np.linspace(18, 29 , num_buckets, retstep=True)
vert_range, vert_step = np.linspace(65, 100, num_buckets, retstep=True)
hori_step = hori_step / 2
vert_step = vert_step / 2
# store each bucket as a tuple in a the buckets dictionary.
buckets = dict()
keys = [(hori, vert) for hori in hori_range for vert in vert_range]
# set each bucket as empty
for key in keys:
    buckets[key] = 0
    
# loop through the randomly sampled data
for month, height in zip(random_months, random_heights):
    # check each bucket and see if randomly sampled data 
    for key in buckets:
        if month > key[0] - hori_step and month < key[0] + hori_step:
            if height > key[1] - vert_step and height < key[1] + vert_step:
                buckets[key] += 1
                break # can only fit in a single bucket
Explanation: Next, we want to get the intensity of the data at locations. We do that by adding the randomly sampled values to buckets. This gives us intensity values for a grid of pixels in our sample range.
End of explanation
pcolor_months = []
pcolor_heights = []
pcolor_densities = []
for key in buckets:
    pcolor_months.append(key[0])
    pcolor_heights.append(key[1])
    pcolor_densities.append(buckets[key])
Explanation: Next, we unpack the buckets into three vectors, pcolor_months, pcolor_heights, and pcolor_densities
End of explanation
def append_to_file(path, data):
    
    append_to_file appends a line of data to specified file.  Then adds new line
    
    Args:
        path (string): the file path
    
    Return:
        VOID
    
    with open(path, 'a') as file:
        file.write(data + '\n')
        
def delete_file_contents(path):
    
    delete_file_contents deletes the contents of a file
    
    Args:
        path: (string): the file path
        
    Return:
        VOID
    
    with open(path, 'w'):
        pass
def threeSequenceCSV(x, y, z):
    
    Writes the x, y, z arrays to a CSV
    
    Args:
        x (sequence): x data
        y (sequence): y data
        z (sequence): z data    
    
    file_name = 'intensityData.csv'
    delete_file_contents(file_name)
    for xi, yi, zi in zip(x, y, z):
        append_to_file(file_name, "{}, {}, {}".format(xi, yi, zi))
        
def twoSequenceCSV(x, y):
    
    Writes the x, y arrays to a CSV
    
    Args:
        x (sequence): x data
        y (sequence): y data
    
    file_name = 'monthsHeights.csv'
    delete_file_contents(file_name)
    
    for xi, yi in zip(x, y):
        append_to_file(file_name, "{}, {}".format(xi, yi))
        
def fittedLineCSV(x, slope, intercept):
    
    Writes line data to a CSV
    
    Args:
        x (sequence): x data
        slope (float): slope of line
        intercept (float): intercept of line
    
    file_name = 'fittedLineCSV.csv'
    delete_file_contents(file_name)
    for xi in x:
        append_to_file(file_name, "{}, {}".format(xi, slope*xi + intercept))
        
def makeCSVData(random_x, random_y, intensities, original_x, original_y, slope, intercept):
    
    Calls the 3 csv making functions with appropriate parameters.
    
    threeSequenceCSV(random_x, random_y, intensities)
    twoSequenceCSV(original_x, original_y)
    fittedLineCSV(original_x, slope, intercept)
makeCSVData(pcolor_months, pcolor_heights, pcolor_densities, ages, heights, slope, intercept)
Explanation: Since density plotting is much simpler in Mathematica, we are going to export all our data to csv files and plot them in Mathematica
End of explanation
def probHeightRange(heights, low, high):
    
    probHeightRange returns the probability that height is within a particular range
    
    Args:
        height (sequence): sequence of heights
        low (number): the bottom of the range
        high (number): the top of the range
        
    Returns:
        prob (float): the probability of being in the height range
    
    successes = 0
    total = len(heights)
    for height in heights:
        if low < height and height < high:
            successes += 1
    
    return successes / total
    
probHeightRange(random_heights, 0, 76.1)
Explanation: Now, we have a figure that contains superimposed density plot, original data scatter plot, and least squares line plot. While the density plot is not strictly necessary, it gives us a much better idea as to how large of a distribution we have. Note that the x axis is ages (months) and the y axis is height (cm).
<img src="ageHeightAllPlots5.png" alt="Density Plot with orignal data/fit" height="400" width="400">
So, now we can go back to our original question. What is the probability that a child between the ages of 18 - 29 is less than 76.1 cm (the lowest height in our dataset)? To answer this question, we sample our random data and see whether or not it is less than 76.1 cm. We will do this using a function.
End of explanation
<END_TASK>