From 081c06fb9b52954f2cd1fa5fe1550ae77ad999d9 Mon Sep 17 00:00:00 2001 From: Patrick Doherty Date: Thu, 23 Nov 2017 23:02:21 -0600 Subject: [PATCH 1/2] Question 2 done --- 12.2.html | 11967 +++++++++++++++++++++++++++++++++++++++++++++++++++ 12.2.ipynb | 159 + 12.2.py | 54 + 3 files changed, 12180 insertions(+) create mode 100755 12.2.html create mode 100755 12.2.ipynb create mode 100755 12.2.py diff --git a/12.2.html b/12.2.html new file mode 100755 index 0000000..d08b3bf --- /dev/null +++ b/12.2.html @@ -0,0 +1,11967 @@ + + + +12.2 + + + + + + + + + + + + + + + + + + + +
+
+ +
+
+
+
+

Using Regular Expressions to search for matching information

    +
  • Regular expressions can be useful to seek out imformation that matches your data set's information
    • +
  • For example, you could type a site name 4 different ways, but want to search for all of those data points instead of searching by each different style of input
    • +
  • Follow along in this exercise to see 3 examples of utilizing Regular Expressions. The first step is to import any necessary packages
    • +
+ +
+
+
+
+
+
In [2]:
+
+
+
import re
+
+ +
+
+
+ +
+
+
+
+
+

Exercise 1: Time

    +
  • Times after noon but before midnight in military time (12:01 - 23:59)
    • +
+ +
+
+
+
+
+
In [3]:
+
+
+
time= "12:01", "09:25","23:30", "11:59"
+reg1=re.compile("(([1][3-9]|[2][0-3]):[0-5][0-9])|(12:[0-5][1-9])")
+print(filter(reg1.match,time))
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
('12:01', '23:30')
+
+
+
+ +
+
+ +
+
+
+
+
+

In this example, we used a regex to filter out any times that do not match our target. (([1][2-9]|[2][0-3]):[0-5][0-9]) matched anything from 13:00 to 23:59 and our third section of (12:[0-5][1-9]) matched anything from 12:01 to 12:59.

+ +
+
+
+
+
+
+
+

Exercise 2: Genus species names

    +
  • Searching for names in the style of H. sapien
    • +
+ +
+
+
+
+
+
In [4]:
+
+
+
names= "M. tuberculosis", "e. Coli", "E. coli", "Staph aureus"
+reg2=re.compile("[A-Z]\\.\\s[a-z]+")
+print(filter(reg2.match,names))
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
('M. tuberculosis', 'E. coli')
+
+
+
+ +
+
+ +
+
+
+
+
+

In this example we used the expression [A-Z]\.\s[a-z]+ to indicate we wanted terms with only 1 capital letter, followed by period, a space, and then one or more undercase letters(indicated by the plus sign at the end of the [a-z]

+ +
+
+
+
+
+
+
+

Exercise 3: Social Security Numbers

    +
  • Searching for numbers in the correct format xxx-xx-xxxx
    • +
+ +
+
+
+
+
+
In [5]:
+
+
+
numbers= "123-12-1234", "Bob", "234-234-234", "999-88-7777", "Using this to find SS numbers to commit fraud is illegal"
+reg3=re.compile("[0-9]{3}-[0-9]{2}-[0-9]{4}")
+print(filter(reg3.match,numbers))
+
+ +
+
+
+ +
+
+ + +
+ +
+ + +
+
('123-12-1234', '999-88-7777')
+
+
+
+ +
+
+ +
+
+
+
+
+

In this last example, we used the expression [0-9]{3}-[0-9]{2}-[0-9]{4} to filter out anything that didn't follow the xxx-xx-xxxx format, including the dashes. Only numbers were allowed to fill the "x" spaces in the example.

+

With these 3 examples, you should have a firmer grasp on the benefits of using Regular Expressions to find terms in a series of strings.

+ +
+
+
+
+
+ + + + + + diff --git a/12.2.ipynb b/12.2.ipynb new file mode 100755 index 0000000..4efb774 --- /dev/null +++ b/12.2.ipynb @@ -0,0 +1,159 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Using Regular Expressions to search for matching information\n", + "* Regular expressions can be useful to seek out imformation that matches your data set's information\n", + "* For example, you could type a site name 4 different ways, but want to search for all of those data points instead of searching by each different style of input\n", + "* Follow along in this exercise to see 3 examples of utilizing Regular Expressions. The first step is to import any necessary packages" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import re" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 1: Time\n", + "* Times after noon but before midnight in military time (12:01 - 23:59) " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "('12:01', '23:30')\n" + ] + } + ], + "source": [ + "time= \"12:01\", \"09:25\",\"23:30\", \"11:59\"\n", + "reg1=re.compile(\"(([1][3-9]|[2][0-3]):[0-5][0-9])|(12:[0-5][1-9])\")\n", + "print(filter(reg1.match,time))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example, we used a regex to filter out any times that do not match our target. (([1][2-9]|[2][0-3]):[0-5][0-9]) matched anything from 13:00 to 23:59 and our third section of (12:[0-5][1-9]) matched anything from 12:01 to 12:59." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercise 2: Genus species names\n", + "* Searching for names in the style of H. sapien" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "('M. tuberculosis', 'E. coli')\n" + ] + } + ], + "source": [ + "names= \"M. tuberculosis\", \"e. Coli\", \"E. coli\", \"Staph aureus\"\n", + "reg2=re.compile(\"[A-Z]\\\\.\\\\s[a-z]+\")\n", + "print(filter(reg2.match,names))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we used the expression [A-Z]\\\\.\\\\s[a-z]+ to indicate we wanted terms with only 1 capital letter, followed by period, a space, and then one or more undercase letters(indicated by the plus sign at the end of the [a-z]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercise 3: Social Security Numbers\n", + "* Searching for numbers in the correct format xxx-xx-xxxx" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "('123-12-1234', '999-88-7777')\n" + ] + } + ], + "source": [ + "numbers= \"123-12-1234\", \"Bob\", \"234-234-234\", \"999-88-7777\", \"Using this to find SS numbers to commit fraud is illegal\"\n", + "reg3=re.compile(\"[0-9]{3}-[0-9]{2}-[0-9]{4}\")\n", + "print(filter(reg3.match,numbers))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this last example, we used the expression [0-9]{3}-[0-9]{2}-[0-9]{4} to filter out anything that didn't follow the xxx-xx-xxxx format, including the dashes. Only numbers were allowed to fill the \"x\" spaces in the example.\n", + "\n", + "With these 3 examples, you should have a firmer grasp on the benefits of using Regular Expressions to find terms in a series of strings. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/12.2.py b/12.2.py new file mode 100755 index 0000000..6ade34b --- /dev/null +++ b/12.2.py @@ -0,0 +1,54 @@ + +# coding: utf-8 + +# # Using Regular Expressions to search for matching information +# * Regular expressions can be useful to seek out imformation that matches your data set's information +# * For example, you could type a site name 4 different ways, but want to search for all of those data points instead of searching by each different style of input +# * Follow along in this exercise to see 3 examples of utilizing Regular Expressions. The first step is to import any necessary packages + +# In[2]: + + +import re + + +# # Exercise 1: Time +# * Times after noon but before midnight in military time (12:01 - 23:59) + +# In[3]: + + +time= "12:01", "09:25","23:30", "11:59" +reg1=re.compile("(([1][3-9]|[2][0-3]):[0-5][0-9])|(12:[0-5][1-9])") +print(filter(reg1.match,time)) + + +# In this example, we used a regex to filter out any times that do not match our target. (([1][2-9]|[2][0-3]):[0-5][0-9]) matched anything from 13:00 to 23:59 and our third section of (12:[0-5][1-9]) matched anything from 12:01 to 12:59. + +# ## Exercise 2: Genus species names +# * Searching for names in the style of H. sapien + +# In[4]: + + +names= "M. tuberculosis", "e. Coli", "E. coli", "Staph aureus" +reg2=re.compile("[A-Z]\\.\\s[a-z]+") +print(filter(reg2.match,names)) + + +# In this example we used the expression [A-Z]\\.\\s[a-z]+ to indicate we wanted terms with only 1 capital letter, followed by period, a space, and then one or more undercase letters(indicated by the plus sign at the end of the [a-z] + +# ## Exercise 3: Social Security Numbers +# * Searching for numbers in the correct format xxx-xx-xxxx + +# In[5]: + + +numbers= "123-12-1234", "Bob", "234-234-234", "999-88-7777", "Using this to find SS numbers to commit fraud is illegal" +reg3=re.compile("[0-9]{3}-[0-9]{2}-[0-9]{4}") +print(filter(reg3.match,numbers)) + + +# In this last example, we used the expression [0-9]{3}-[0-9]{2}-[0-9]{4} to filter out anything that didn't follow the xxx-xx-xxxx format, including the dashes. Only numbers were allowed to fill the "x" spaces in the example. +# +# With these 3 examples, you should have a firmer grasp on the benefits of using Regular Expressions to find terms in a series of strings. From aa8f8e4d1e118bda34e77de6181b9731bc144b28 Mon Sep 17 00:00:00 2001 From: Tim Burton Date: Fri, 24 Nov 2017 01:30:27 -0500 Subject: [PATCH 2/2] exercise 1 done --- 12.1.html | 12555 +++++++++++++++++++++++++++++++++++++++++++++++++++ 12.1.ipynb | 530 +++ 12.1.py | 108 + 3 files changed, 13193 insertions(+) create mode 100644 12.1.html create mode 100644 12.1.ipynb create mode 100644 12.1.py diff --git a/12.1.html b/12.1.html new file mode 100644 index 0000000..8c3bc41 --- /dev/null +++ b/12.1.html @@ -0,0 +1,12555 @@ + + + +12.1 + + + + + + + + + + + + + + + + + + + +
+
+ +
+
+
+
+
+

General description of "chickwts.txt" dataset.

Import packages and data

+
+
+
+
+
+
In [20]:
+
+
+
import pandas as pd
+import numpy as np
+from scipy.stats import chi2
+from scipy.optimize import minimize
+from scipy.stats import norm
+from plotnine import *
+
+chickens=pd.read_table("chickwts.txt", delimiter=",")
+
+ +
+
+
+ +
+
+
+
+
+
+

1. Calculate mean weight of chickens based on feed source

create an empty dictionary for "means". Fill by using a for loop to loop through each type of chicken feed, and add the mean weight of chickens for each feed to the dictionary entry for that feed.

+
+
+
+
+
+
In [21]:
+
+
+
means=dict()
+
+for i in np.unique(chickens.feed):
+    means[i]=np.mean(chickens[chickens.feed == i][["weight"]])[0]
+    
+means_df=pd.DataFrame(means, index =[0])
+means_df
+
+ +
+
+
+ +
+
+ + +
+
Out[21]:
+ + +
+
+ + + + + + + + + + + + + + + + + + + + + + + + +
caseinhorsebeanlinseedmeatmealsoybeansunflower
0323.583333160.2218.75276.909091246.428571328.916667
+
+
+ +
+ +
+
+ +
+
+
+
+
+
+

2. Generate a plot summarizing chicken weight data

use ggplot2 to create a boxplot of "chickwts.txt" data, plotting chicken weight on the y axis, grouped by chicken feed

+
+
+
+
+
+
In [22]:
+
+
+
ggplot(chickens, aes(x="feed",y="weight")) + geom_boxplot() + theme_classic()
+
+ +
+
+
+ +
+
+ + +
+
+ + + +
+ +
+ +
+ +
+
Out[22]:
+ + + +
+
<ggplot: (8741302840413)>
+
+ +
+ +
+
+ +
+
+
+
+
+
+

3. Null and alternative hyptheses for the difference in chick weight when fed soybean feed vs sunflower feed

Null hypothesis: mean(sunflower)-mean(soybean) = 0

Alternative hypothesis: mean(sunflower)-mean(soybean) =/= 0

+
+
+
+
+
+
+
+
+

4. Test Null and alternative hypotheses using chi-2 likelihood test

Create a data "sun_soy_chickens" containing only data of interest (chickens fed on soybean or sunflower feed). Add a column "Factor" to this new data frame, containing a 1 for chickens fed with sunflower feed, and a 0 for chickens fed with soybean feed.

+
+
+
+
+
+
In [23]:
+
+
+
sun_soy_chickens=chickens[(chickens.feed == "sunflower") | (chickens.feed == "soybean")]
+
+sun_soy_chickens["Factor"]=0
+for i in range(0, len(sun_soy_chickens.feed)):
+    if sun_soy_chickens.feed.iloc[i] == "sunflower":
+        sun_soy_chickens.iloc[i,2]=1
+    else:
+        sun_soy_chickens.iloc[i,2]=0
+        
+sun_soy_chickens
+
+ +
+
+
+ +
+
+ + +
+
+ +
+
/home/tb/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
+A value is trying to be set on a copy of a slice from a DataFrame.
+Try using .loc[row_indexer,col_indexer] = value instead
+
+See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
+  This is separate from the ipykernel package so we can avoid doing imports until
+
+
+
+ +
+
Out[23]:
+ + +
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
weightfeedFactor
22243soybean0
23230soybean0
24248soybean0
25327soybean0
26329soybean0
27250soybean0
28193soybean0
29271soybean0
30316soybean0
31267soybean0
32199soybean0
33171soybean0
34158soybean0
35248soybean0
36423sunflower1
37340sunflower1
38392sunflower1
39339sunflower1
40341sunflower1
41226sunflower1
42320sunflower1
43295sunflower1
44334sunflower1
45322sunflower1
46297sunflower1
47318sunflower1
+
+
+ +
+ +
+
+ +
+
+
+
+
+
+

Define null and alternative hypothesis functions

+
+
+
+
+
+
In [24]:
+
+
+
def nllalt(p,obs):
+    B0=p[0]
+    sigma=p[1]
+    B1=p[2]
+    
+    expected=B0+B1*obs.Factor
+    nll=-1*norm(expected,sigma).logpdf(obs.weight).sum()
+    return nll
+    
+    
+def nllnull(p,obs):
+    B0=p[0]
+    sigma=p[1]
+    
+    expected=B0
+    nll=-1*norm(expected,sigma).logpdf(obs.weight).sum()
+    return nll
+
+ +
+
+
+ +
+
+
+
+
+
+

Make guess, use minimize function to find optimal coefficients for alternative and null hypothesis functions. T test to determine the significance of the factor (presence of sunflower seed) in determining the weight of chickens in this study

+
+
+
+
+
+
In [25]:
+
+
+
guess=[1,1,1]
+
+fit_alt=minimize(nllalt, guess, method="Nelder-Mead", options={'disp': True},args=sun_soy_chickens)
+fit_null=minimize(nllnull, guess, method="Nelder-Mead", options={'disp': True},args=sun_soy_chickens)
+
+D=2*(fit_null.fun-fit_alt.fun)
+1-chi2.cdf(x=D, df=1)
+
+ +
+
+
+ +
+
+ + +
+
+ +
+
Optimization terminated successfully.
+         Current function value: 138.469162
+         Iterations: 200
+         Function evaluations: 363
+Optimization terminated successfully.
+         Current function value: 145.240592
+         Iterations: 138
+         Function evaluations: 261
+
+
+
+ +
+
Out[25]:
+ + + +
+
0.00023317672869682671
+
+ +
+ +
+
+ +
+
+
+
+
+
+

Interpreting the result of the likelihood ratio test

Since the P-value is \<0.05, there is a significant difference in chicken weights between chickens fed on sunflower and soybean feed.

+
+
+
+
+
+ + + + + + diff --git a/12.1.ipynb b/12.1.ipynb new file mode 100644 index 0000000..d5e585f --- /dev/null +++ b/12.1.ipynb @@ -0,0 +1,530 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# General description of \"chickwts.txt\" dataset.\n", + "\n", + "\n", + "\n", + "### Import packages and data" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "from scipy.stats import chi2\n", + "from scipy.optimize import minimize\n", + "from scipy.stats import norm\n", + "from plotnine import *\n", + "\n", + "chickens=pd.read_table(\"chickwts.txt\", delimiter=\",\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1. Calculate mean weight of chickens based on feed source\n", + "\n", + "#### create an empty dictionary for \"means\". Fill by using a for loop to loop through each type of chicken feed, and add the mean weight of chickens for each feed to the dictionary entry for that feed." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
caseinhorsebeanlinseedmeatmealsoybeansunflower
0323.583333160.2218.75276.909091246.428571328.916667
\n", + "
" + ], + "text/plain": [ + " casein horsebean linseed meatmeal soybean sunflower\n", + "0 323.583333 160.2 218.75 276.909091 246.428571 328.916667" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "means=dict()\n", + "\n", + "for i in np.unique(chickens.feed):\n", + " means[i]=np.mean(chickens[chickens.feed == i][[\"weight\"]])[0]\n", + " \n", + "means_df=pd.DataFrame(means, index =[0])\n", + "means_df" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2. Generate a plot summarizing chicken weight data\n", + "\n", + "#### use ggplot2 to create a boxplot of \"chickwts.txt\" data, plotting chicken **weight** on the y axis, grouped by chicken **feed**" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAFzCAYAAAB4qqApAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl8VNX9//H3TPaFGAJMEkUSZA8gEioqi4AgtVIQhVK2\nypfI7veB+FOxqIDUKn0EFSigRgRB2awVt1o1VTZBESWg7LKGCCEBQoAsk2Qy5/eHX6ZEQCDOnQnh\n9Xw8eOjcOzPnc24md9455y42Y4wRAACABez+LgAAAFRfBA0AAGAZggYAALAMQQMAAFiGoAEAACxD\n0AAAAJYhaAAAAMsQNAAAgGUIGgAAwDIEDQAAYJmrLmgUFRUpIyNDRUVF/i4FAIBq76oLGjt37lSb\nNm20c+dOf5cCAEC1d9UFDQAA4DsEDQAAYBmCBgAAsAxBAwAAWIagAQAALEPQAAAAliFoAAAAyxA0\nAACAZQgaAADAMgQNAAD8wO12q6SkxN9lWI6gAQCADxlj9Morr6hdu3Zq3769+vXrp8zMTH+XZRmC\nBgAAPrR06VK9/vrrcrlckqTMzEyNHj262t7sk6ABAIAPffjhhyovL/c8Li8v19GjR7V9+3Y/VmUd\nggYAAD5kjLms5Vc6ggYAAD509913KyAgwPPYbrerVq1aSkpK8mNV1iFoAADgQ4MHD9agQYNkt//0\nFVy3bl29/PLLioiI8HNl1rCZ6jpWcwEZGRlq06aNNm7cqOTkZH+XAwC4SrlcLjmdTkVGRvq7FEsx\nogEAgB8EBgZW+5AhETQAAICFCBoAAMAyBA0AAGAZggYAALAMQQMAAFiGoAEAACxD0AAAAJYhaAAA\nAMsQNAAAgGUIGgAAwDIEDQAAYJlAfxdwtlOnTmn06NGKj4/X888/L0nKzMzUrFmzdODAAcXGxmrE\niBFq1aqV5zXr1q3TwoULlZeXp6ZNm2rs2LFyOBz+6gIAADhLlRrReP3113X99dd7HrtcLj3zzDNq\n27atli5dqv79+2vq1KnKz8+XJGVlZWnmzJkaPXq0Fi9erMTERKWmpvqrfMsUFBQoLS1NBQUF/i4F\nAIDLUmWCxtatW3X48GF169bNs2zLli0qKSlR3759FRQUpI4dOyohIUHr1q2TJK1atUrJyclq3bq1\nQkJCNHDgQB04cEAHDx70VzcsUVBQoLlz5xI0AABXnCoRNMrKypSWlqZRo0bJZrN5lh88eFCJiYmy\n2/9bZv369ZWZmSnpp2mV+vXre9aFh4crLi7Os/6M7OxsZWRkKCMjQzt27LC4NwAA4IwqcYzGO++8\no1atWql+/frat2+fZ3lxcbEiIiIqPDciIkK5ubmSJKfTec768PBwFRcXV1iWlpamKVOmWFQ9AAC4\nEL+PaBw+fFiff/65Bg4ceM66sLAwFRYWVlhWWFiosLAwSVJoaKiKioouuP6MkSNHauPGjdq4caMW\nLVrk5R4AAIAL8fuIxo4dO3TixAmNGjVKklRaWqrS0lLdf//9GjNmjDIzM+V2uz3TJ/v27VOnTp0k\nSQkJCRVGQIqKipSTk6OEhIQKbcTHxys+Pt5HPQIAAGf4PWh06NBBycnJnsdffPGFVq5cqUmTJqlG\njRoKDg7W8uXLdc899+jrr7/WwYMH1b59e0lS586d9cgjj2jz5s1KSkrSkiVLlJiYqHr16vmrOwAA\n4Cx+DxohISEKCQnxPI6IiFBAQIBq1qwpSXrqqac0e/ZsLVu2TA6HQxMmTFB0dLQk6frrr9fYsWM1\nZ84cnThxQk2aNNH48eP90g8AAC7VyZMn9Z///EenT5/WjTfeqDZt2vi7JMv4PWj8XNeuXdW1a1fP\n48TERM/Fu86nQ4cO6tChgy9KAwDgVzty5IiGDBmiU6dOyWazqaysTP/7v/+rIUOG+Ls0S/j9YFAA\nAK4mL7zwgk6ePKmysjKVlpbKGKNZs2YpKyvL36VZgqABAIAP7d27Vy6Xq8Iyu91e7S42eQZBAwAA\nH4qPj69wIUpJcrvdio2N9VNF1iJoAADgQw8//LCCgoIUGBgom82mgIAA3XvvvWrYsKG/S7NElTsY\nFACA6qxhw4ZaunSpli9frtOnT6tVq1bq2bOnv8uyDEEDAAAfq1evnsaNG+fvMnyCqRMAAGAZggYA\nALAMQQMAAFiGoAEAACxD0AAAAJYhaAAAAMsQNADgIgoKCpSWlqaCggJ/lwJccQgaAHARBQUFmjt3\nLkEDqASCBgAAsAxBAwAAWIagAQAALEPQAAAAliFoAAAAy3D31krKz8+X0+n0SVu5ubkV/usLoaGh\nio6O9ll7AIDqiaBRCfn5+erevbvcbrdP201JSfFZW3a7Xenp6YQNAMCvQtCoBKfTKbfbrfnz58vh\ncPi7HK/Lzc1VSkqKz0ZsAADVF0HjV3A4HIqLi/N3GQAAVFkcDAoAACxD0AAAAJYhaAAAAMsQNAAA\ngGUIGgAAwDIEDQAAYBmCBgAAsAxBAwAAWIagAQAALEPQAAAAliFoAAAAyxA0AACAZQgaAADAMgQN\nAABgGYIGAACwDEEDAABYhqABAAAsQ9AAAACWIWgAAADLEDQAAIBlCBoAAMAygf4u4Eq2fft2HTly\nxN9leF1eXp6/SwAAVBMEjV9h/Pjx/i4BAIAqjakTAABgGUY0foXU1FTFxMT4uwyvy8vLY7QGwFUp\nPz9fTqfTJ20VFhbqvffeU+/evRUREeGTNkNDQxUdHe2Tts4gaPwKSUlJiouL83cZXlcdjzsBgIvJ\nz89X9+7d5Xa7fdru0qVLfdaW3W5Xenq6T8MGQQMAAElOp1Nut1vz58+Xw+Hwdzlel5ubq5SUFJ+N\n2JxB0AAA4CwOh6Najlb7CweDAgAAyxA0AACAZQgaAADAMgQNAABgGYIGAACwDGedAABwFu5j5V0E\nDQAAzsKVkb2LqRMAgNcUFBQoLS1NBQUF/i4FVQQjGgAArykoKNDcuXN1zz33KDIy0t/lVAr3sfIu\nggYAAGfhPlbexdQJAACwDEEDAABYhqABAAAsQ9AAAACW4WBQAFek/Px8OZ1On7SVm5tb4b++EBoa\nqujoaJ+1B1iFoAHgipOfn6/u3bvL7Xb7tN2UlBSftWW325Wenk7YwBWPoAHgiuN0OuV2uzV//nw5\nHA5/l+N1ubm5SklJ8dmIDWAlggaAK5bD4aiW1zsAqhMOBgUAAJYhaAAAAMswdQIAwFl8eXaRL/mr\nXwQNAAD00ynFdrvdp2cX+ZrdbldoaKhP2yRoAAAgKTo6Wunp6T69PktKSopPz57yx/VZCBoAAPwf\nf1y3pLqfPcXBoMBVqqCgQGlpaSooKPB3KQCqMYIGcJUqKCjQ3LlzCRoALEXQAAAAliFoAAAAy1SJ\ng0Fnz56tb7/9VsXFxapRo4a6d++ufv36SZIyMzM1a9YsHThwQLGxsRoxYoRatWrlee26deu0cOFC\n5eXlqWnTpho7dmy1vPcBAABXoioxotGrVy+9/PLLeuuttzR16lStXr1aa9eulcvl0jPPPKO2bdtq\n6dKl6t+/v6ZOnar8/HxJUlZWlmbOnKnRo0dr8eLFSkxMVGpqqp97AwAAzqgSQaNevXoKCwvzPLbZ\nbDp8+LC2bNmikpIS9e3bV0FBQerYsaMSEhK0bt06SdKqVauUnJys1q1bKyQkRAMHDtSBAwd08OBB\nf3UFAACcpUpMnUjSwoUL9a9//UslJSVyOBzq0qWLvvzySyUmJspu/28eql+/vjIzMyX9NK3SqFEj\nz7rw8HDFxcUpMzNT9erV8yzPzs5Wdna2JGnHjh0+6hEAAKgyQWPIkCG6//77tWfPHq1fv14REREq\nLi5WREREhedFRER4rtfudDrPWR8eHq7i4uIKy9LS0jRlyhRrOwAAAM5RJaZOzrDZbGrUqJGCgoK0\ndOlShYWFqbCwsMJzCgsLPdMsoaGhKioquuD6M0aOHKmNGzdq48aNWrRokbWdAAAAHlVmRONsbrdb\n2dnZSk5O1jvvvCO32+2ZPtm3b586deokSUpISNC+ffs8rysqKlJOTo4SEhIqvF98fLzi4+N91wEA\nACCpCoxoFBQUaOXKlSoqKpLb7db27dv18ccf66abblLLli0VHBys5cuXq6ysTGvXrtXBgwfVvn17\nSVLnzp2VkZGhzZs3q7S0VEuWLFFiYmKF4zMAAID/VIkRjc8++0yvvvqq3G63YmJi1Lt3b/Xo0UM2\nm01PPfWUZs+erWXLlsnhcGjChAmem95cf/31Gjt2rObMmaMTJ06oSZMmGj9+vJ97AwAAzvB70IiM\njNSzzz57wfWJiYl6/vnnL7i+Q4cO6tChgxWlAQBgmcjISA0fPlyRkZH+LsVSlZo6WbNmzQVvxFRQ\nUKA1a9b8qqIAAKjuIiMjNXLkSILG+XTp0kXbt28/77pdu3apS5cuv6ooAABQPVQqaBhjLrjufKeX\nAgCAq9MlH6Oxfv16ffnll57HS5Ys0dq1ays8x+l06v3331ezZs28VyEAALhiXXLQ+PTTTz1X17TZ\nbPr73/9+znOCgoLUrFkzvfTSS96rEAAAXLEueepk8uTJcrvdcrvdMsZo/fr1nsdn/pWUlGjz5s1q\n166dlTUDAIArRKVOb3W73d6uAwAAVEOVvo5GeXm5vv76a/34449yOp3nrL///vt/VWEAAODKV6mg\nkZGRofvuu09ZWVnnPQPFZrMRNAAAQOWCxujRo3XNNddo4cKFSkpKUnBwsLfrAgAA1UClgsa2bdv0\n9ttve+6iCgAAcD6VumBX48aNderUKW/XAgAAqplKBY3p06dr6tSp2rlzp7frAQAA1cglT520bNlS\nNpvN8zg7O1stWrTQtdde67lt+xk2m03fffed96oEAABXpEsOGm3atKkQNCDl5ub6uwRLVNd+AQB8\n75KDxoIFCyws48oSGhoqu92ulJQUf5diGbvdrtDQUH+XAQC4wlX6gl1Xs+joaKWnp5/3QmVWyM3N\nVUpKiubPny+Hw+GTNkNDQ8+ZEgMA4HJVKmj80l/ydrtd11xzjVq3bq377rtP4eHhlS6uKvPHl7DD\n4VBcXJzP2wUAoLIqFTQ2bdqkw4cP6+jRo4qJiZHD4VBubq7y8vJUp04dRUREaObMmXryySe1YsUK\nNWjQwNt1AwCAK0ClTm+dNm2aoqKi9MUXX+jYsWPavn27jh07ptWrVysqKkpz5szRjh07FBISovHj\nx3u7ZgAAcIWo1IjGo48+qqefflrt27evsLxjx46aNGmSHnvsMW3dulUTJkzQI4884pVCAQCVk5+f\n79Njys7+ry9wTFnVVqmgsWvXrgv+UGvWrKm9e/dKkho0aKDi4uLKVwcA+FXy8/PVvXt3ud1un7br\ny7Py7Ha70tPTCRtVVKWCRtOmTfX888+rS5cuFQ72LCws1LRp05SUlCRJOnz4sGJjY71TKQDgsjmd\nTrndbp+eteZLZ87K89WIDS5fpYLGrFmz9Lvf/U5169ZVly5dVKdOHR09elQrVqyQy+XSJ598Ikn6\n/vvv1bdvX68WDAC4fJy1Bn+pVNDo0KGDdu/erRdffFHffvuttm/frvj4eI0YMUIPP/yw58P83HPP\nebVYAABwZan0Bbvi4uKUmprqzVoA4LJs375dR44c8XcZXpeXl+fvEgCv4cqgAK5YnD4PVH2XHDRu\nvPFGLVmyRC1atDjnTq4/x91bAQCAdJl3b42IiPD8P3dyBeBvqampiomJ8XcZXpeXl8doDaqNSw4a\nr7/+uuf/uZMrgKogKSmpWp5JUR2PO8HVq1KXID+bMUaHDx+Wy+XyRj0AAKAaqfTBoJ9++qkmT56s\nTZs2yeVy6ZtvvlFycrJGjBihTp06adCgQd6sE7gqcKloANVNpYLG0qVLNXjwYPXr10/Dhw/X8OHD\nPesaNGig119/naABXCYuFQ2gOqpU0HjmmWc0btw4vfDCCyovL68QNJo3b67p06d7rUDgasGlogFU\nR5UKGvv27dPdd9993nURERE6efLkryoKuJpxqWgA1UmlDgaNi4vTzp07z7vu+++/V0JCwq8qCgAA\nVA+VChoDBw7U008/rc8//9yzzGazaevWrUpNTdXgwYO9ViAAALhyVWrq5Omnn9a2bdt05513qlat\nWpKk3/3udzp69Kh+//vf689//rNXiwQAAFemSgWN4OBgvf/++1q5cqXS09N1/PhxxcTEqFu3burW\nrZu3awSA8/Llqbm+VF37hatTpYJGz5491alTJ3Xs2FF//etfFRAQ4O26AOCCQkNDZbfbfXpqrq/Z\n7XaFhoZ67f240y38pVJBIzIyUjNmzND48eMVERGhW2+9VR07dtTtt9+u2267TSEhId6uEwA8oqOj\nlZ6e7tOLm6WkpPj01GNvX9yMe6fAXyp9wS5J2rNnj9asWaMvvvhCCxcu1JQpUxQUFKSbb75ZX3zx\nhVcLBYCz+eOiX5x6DFy+Sl+CXJIaNmyohg0bqlOnTlq1apUWL16sVatW6csvv/RWfQAAL+BOt/CX\nSgWNHTt2aM2aNVq9erXWrFmjI0eOqHnz5rr99ts1evRo3X777d6uEwDwK3CnW/hLpYJG8+bNFRYW\npiFDhujll19Whw4dVLNmTW/XBgAArnCVCho9evTQunXrNG/ePH3//ff66quv1KlTJ7Vv316RkZHe\nrhEAAFyhKnVl0A8//FDHjx/Xhg0b9Mc//lG7d+/WkCFDFBMTo5tvvlmPPvqot+sEAABXoEofDGqz\n2dSqVSu1atVKffv21erVq5WWlqbVq1crIyNDzz//vDfrBAAAV6BKBY39+/drzZo1nn/79u1TUFCQ\nWrdurccee0ydOnXydp0AAOAKVKmg0aBBA4WGhqpt27YaMGCAOnXqpNtuu03h4eHerg8AAFzBKhU0\n1qxZo7Zt2yo4ONjb9QAAgGqkUkGjQ4cO3q4DAABUQ5U66wQAAOBSEDQAAIBlCBoAAMAyBA0AAGAZ\nggYAALAMQQMAAB9777331KNHD3Xu3Fljx47V8ePH/V2SZQgaAAD40EcffaTnnntOOTk5Kigo0IYN\nGzRq1CiVlpb6uzRLEDRQrRQUFCgtLU0FBQX+LgUAzmvZsmVyu92exy6XS/v379f27dv9WJV1CBqo\nVgoKCjR37lyCBoAqq6Sk5JxlNpvtvMurA4IGAAA+dPvttysw8L8X5rbZbAoPD1fTpk39WJV1CBoA\nAPjQyJEjdccdd3ge16hRQzNnztQ111zjx6qsU6l7nQCwzvbt23XkyBF/l+F1eXl5/i4BqBKCgoL0\n3HPPady4cSooKNB1112nkJAQf5dlGYIGUMWMHz/e3yUA8AGHwyGHw+HvMizH1AkAALAMIxpAFZOa\nmqqYmBh/l+F1eXl5jNYAVyGCBlDFJCUlKS4uzt9leF11PO4EwMUxdQIAACxD0LgCREZGavjw4YqM\njPR3KQAAXBamTq4AkZGRGjlypL/LAADgsjGiAQAALEPQAAAAliFoAAAAyxA0AACAZQgaAADAMpx1\nAgBXgdzcXH+XYInq2q/qhKABANVYaGio7Ha7UlJS/F2KZex2u0JDQ/1dBi7A70GjrKxMr7zyir77\n7judPn1atWvXVr9+/dSpUydJUmZmpmbNmqUDBw4oNjZWI0aMUKtWrTyvX7dunRYuXKi8vDw1bdpU\nY8eOvSruhgcAlyI6Olrp6elyOp0+aS83N1cpKSmaP3++z/bFoaGhio6O9klbuHx+Dxrl5eWKiYnR\nX//6V8XGxmrHjh36y1/+otjYWDVs2FDPPPOMunfvrqlTp2r9+vWaOnWqXnnlFUVHRysrK0szZ87U\nhAkTlJSUpDfffFOpqal6/vnn/d0tAKgy/PEl7HA4quU9e3D5/H4waGhoqAYNGqS4uDjZbDYlJSWp\nWbNm2rFjh7Zs2aKSkhL17dtXQUFB6tixoxISErRu3TpJ0qpVq5ScnKzWrVsrJCREAwcO1IEDB3Tw\n4EE/9woAAEhVIGj8nNPp1J49e5SQkKCDBw8qMTFRdvt/y6xfv74yMzMl/TStUr9+fc+68PBwxcXF\nedYDAAD/8vvUydncbrdmzJihRo0aqXXr1vrhhx8UERFR4TkRERGeo4ydTuc568PDw1VcXFxhWXZ2\ntrKzsyVJO3bssLAHAADgbFUmaBhj9NJLLykvL09TpkyRzWZTWFiYCgsLKzyvsLBQYWFhkn6adikq\nKrrg+jPS0tI0ZcoUazsAAADOUSWmTowxeuWVV7R//349/fTTnqBQr149ZWZmyu12e567b98+JSQk\nSJISEhK0b98+z7qioiLl5OR41p8xcuRIbdy4URs3btSiRYt80CMAACBVkaCRlpamXbt2acqUKQoP\nD/csb9mypYKDg7V8+XKVlZVp7dq1OnjwoNq3by9J6ty5szIyMrR582aVlpZqyZIlSkxMVL169Sq8\nf3x8vJKTk5WcnKxmzZr5tG8AAFzN/D51kpubq3//+98KCgqqcEGZvn37ql+/fnrqqac0e/ZsLVu2\nTA6HQxMmTPCcqnX99ddr7NixmjNnjk6cOKEmTZpo/Pjx/uoKAAD4Gb8HDYfDoQ8++OCC6xMTE3/x\nuhgdOnRQhw4drCgNAAD8SlVi6gQAAFRPBA0AAGAZggYAALAMQQMAAFiGoAEAACxD0AAAAJYhaAAA\nAMsQNAAAgGUIGgAAwDIEDQAAYBmCBgAAsAxBAwAAWIagAQAALEPQAAAAlvH7beIBVJSbm+vvEixR\nXfsF4JcRNGC5/Px8OZ1On7R15svMl19qoaGhio6O9sr72O12paSkeKGqqslutys0NNTfZQDwIYIG\nLJWfn6/u3bvL7Xb7tF1fflnb7Xalp6f/6rARHR2t9PR0n4aylJQUzZ8/Xw6HwydteiuUAbhyEDRg\nKafTKbfb7dMvM18682XtrXDgjy9hh8OhuLg4n7cL4OpA0IBP8GUGAFcnzjoBAACWIWgAAADLEDQA\nAIBlCBoAAMAyBA0AAGAZggYAALAMQQMAAFiGoAEAACxD0AAAAJYhaAAAAMsQNAAAgGUIGgAAwDIE\nDQAAYBmCBgAAsAxBAwDgNZGRkRo+fLgiIyP9XQqqiEB/FwAAqD4iIyM1cuRIf5eBKoQRDQAAYBlG\nNOAT27dv15EjR/xdhtfl5eX5uwQAqNIIGvCJ8ePH+7sEoNI47gCoPIIGAFwExx0AlUfQgE+kpqYq\nJibG32V4XV5eHqM1APALCBrwiaSkJMXFxfm7DK+rjsedAIA3cdYJAACwDEEDAABYhqABAAAsQ9AA\nAACWIWgAAADLEDQAAIBlCBoAAMAyBA3gKsVltQH4AhfsAq5SXFYbgC8QNOATubm5/i7BEtW1XwDg\nLQQNWCo0NFR2u10pKSn+LsUydrtdoaGh/i4DAKokggYsFR0drfT0dDmdTp+0l5ubq5SUFM2fP18O\nh8MnbYaGhio6OtonbQHAlYagAcv540vY4XBUy5u4AcCVhrNOAACAZQgaAADAMgQNAABgGYIGAACw\nDEEDAABYhqABAAAsQ9AAAACWIWgAAADLEDQAAIBlCBoAAMAyBA0AAGAZggYAALAMQQMAAFiGoAEA\nACxD0AAAAJYhaAAAAMsQNAAAgGUIGgAAwDIEDQAAYBmCBgAAsAxBAwAAWIagAQAALEPQAAAAliFo\noFqJjIzU8OHDFRkZ6e9SAACSAv1dAOBNkZGRGjlypL/LAAD8H78HjX/9619asWKFDhw4oNtuu02P\nPfaYZ11mZqZmzZqlAwcOKDY2ViNGjFCrVq0869etW6eFCxcqLy9PTZs21dixY+VwOPzRDQAAcB5+\nnzqJiYlRv3791L179wrLXS6XnnnmGbVt21ZLly5V//79NXXqVOXn50uSsrKyNHPmTI0ePVqLFy9W\nYmKiUlNT/dEFAABwAX4PGu3atdOtt96qqKioCsu3bNmikpIS9e3bV0FBQerYsaMSEhK0bt06SdKq\nVauUnJys1q1bKyQkRAMHDtSBAwd08OBBf3QDAACch9+DxoUcPHhQiYmJstv/W2L9+vWVmZkp6adp\nlfr163vWhYeHKy4uzrMeAAD4n9+P0biQ4uJiRUREVFgWERGh3NxcSZLT6TxnfXh4uIqLi895r+zs\nbGVnZ0uSduzYYVHFAADg56ps0AgLC1NhYWGFZYWFhQoLC5MkhYaGqqio6ILrz5aWlqYpU6ZYVywA\nADivKjt1Uq9ePWVmZsrtdnuW7du3TwkJCZKkhIQE7du3z7OuqKhIOTk5nvVnGzlypDZu3KiNGzdq\n0aJF1hcPAAAkVYGgUV5ertLSUrndbrndbpWWlsrlcqlly5YKDg7W8uXLVVZWprVr1+rgwYNq3769\nJKlz587KyMjQ5s2bVVpaqiVLligxMVH16tU7p434+HglJycrOTlZzZo183UXAQC4atmMMcafBSxZ\nskTLli2rsOyOO+7QuHHjdODAAc2ePVsHDhyQw+HQyJEjK1xHY+3atVq4cKFOnDihJk2a6KGHHrro\ndTQyMjLUpk0bbdy4UcnJyZb0CQAA/MTvQcPXCBoAAPiO36dOAABA9VVlzzqxypnTXznNFQCAymna\ntKnCw8Mv6blXXdA4cOCAJGnw4MH+LQQAgCvU5Rx+cNUdo3Hs2DF9+umnSkxMPO81N6qiHTt2aPDg\nwVq0aBFnzVwE2+rSsa0uHdvq0rGtLt2VvK0Y0fgFtWvX1qBBg/xdRqU0a9aMA1gvEdvq0rGtLh3b\n6tKxrS5ddd9WHAwKAAAsE/D0008/7e8icHGRkZHq3LmzatSo4e9Sqjy21aVjW106ttWlY1tduqth\nW111x2gAAADfYeoEAABYhqABAAAsQ9CoRl566SUtWbLE32VIkoYNG6aNGzf6tYYlS5Zo2rRpfq3B\nm85s03/84x+aMWOGv8upYOHChVWupitdr1699OOPP/q7jEsyY8YMLVy40N9lVFknT57Uk08+qT/+\n8Y+aNWtWtds3XcxVd3prdTZmzBh/lwAf6Nevn79LuCoMGzZMo0ePVps2bfxdCq5wn376qUJDQ7Vs\n2TLZbLYq8wehrzCigSqtvLzc3yUAwK+Sk5OjevXqyWaz+bsUSb7frzKiYbHjx49r/vz5+v777+Vy\nudSiRQuvKrynAAAWLklEQVQ9+eSTmjZtmrZu3aqSkhIlJiZq1KhRSkxMlPTTpV0XLFig3NxchYSE\nqEuXLho6dKgkaffu3Zo3b54yMzNVs2ZNDR48WO3atZP00/BlzZo1NWTIEG3ZskXTpk1Tv3799Pbb\nb8vtduvee+/Vfffd57O+Hzx4UIsWLdLhw4eVlJSkRx55RJGRkdq4caMWLlyo3Nxc1a1bV8OGDVPT\npk09fQgODlZ+fr42b96sBx98UNdee63S0tKUlZWloKAgtW7dWo888ogk6fDhw3r11Ve1e/duRURE\nqHfv3rr77rs9NZSVlemFF17Qhg0bVLt2bY0ePVotWrSQJBUVFWnBggX65ptvVF5ernbt2iklJUXB\nwcEqKirSCy+8oB9++EEul0tNmzbVmDFjVKdOHUnSE088oaSkJO3YsUN79uxRQkKCHn30UTkcDsu3\n65IlS3To0CE99thjysnJ0fDhw/Xwww9r8eLFKiwsVNeuXTVs2DBJ0pEjRzRr1izt3btXAQEBuv76\n6/W3v/1NknTixAnNnTtXW7ZsUVBQkLp166b+/fvLbv/p74+VK1fqnXfe0fHjx5WQkKAxY8aoXr16\nkqT9+/dr1qxZOnTokJo3b65atWpZ3u8zhg0bprvvvltr1qzRoUOH1Lp1a40dO1Zz587V+vXr5XA4\n9OijjyohIeEX+3jkyBHNnj1b+/fvlyTddNNNGj16tCIjIzVt2jQdPXpUU6dOld1uV8+ePfWnP/1J\nvXr10ujRo/XBBx/o2LFj6tKliwYNGqSZM2dq27ZtSkhI0Pjx41W7dm1Jv/z53L17t+bOnev5XN92\n220aNmyYgoKCLN+G7777rj788EMVFhYqKipKf/rTn9SxY0ctX75cn3zyiYqKipSUlKRRo0apVq1a\nevfdd/X9999r8uTJnvd45513tG3bNk2aNEmSVFBQoClTpmj79u2qW7euxo4dq4SEBEm//Fn7pZ/D\nmZ93jx49tGbNmnP2Jb7aNj/++KPnd06SSktL1bdvX82dO1exsbGaMWOGQkJCdOLECX333XdyOBx6\n+OGHdcMNN+jFF1/UF198IZvNpo8++kjjxo07p80L7ROPHTum0aNHa8mSJQoKCtKiRYv0zjvvaMmS\nJQoLC9M///lP/fjjjxo3bpzKysq0dOlSrVmzRsXFxWrdurVGjRqlyMhIz37ioYce0tKlSxUSEqI5\nc+ZYsv3Oy8AyLpfLjBs3zrz00kumsLDQlJWVmS1bthhjjPnPf/5jCgsLTWlpqXnttdfMgw8+6Hnd\n/fffb1asWGGMMaaoqMjs3LnTGGPM8ePHzcCBA81XX31lXC6X2blzpxkwYIA5ePCgMcaY6dOnmwUL\nFhhjjPn+++/NPffcY+bPn29KS0vNrl27zL333msOHTrkk74/8MADZuzYsSY3N9cUFRWZxx57zCxa\ntMgcOnTI9OnTx2zYsMG4XC7z+eefm/79+5uTJ096+vCHP/zBfPfdd8btdhun02keffRR89Zbb5ny\n8nJTUlJitm3bZowxxul0mqFDh5qPPvrIlJWVmaysLDN06FCTkZFhjDFm8eLF5p577jErVqwwLpfL\nfPbZZ6Z///7m9OnTxhhjnnvuOTN9+nRTWFhoTp8+bSZNmmTeeOMNY4wxp0+fNmvXrjVOp9MUFRWZ\n1NRUM2XKFE//JkyYYFJSUsyBAwdMaWmpmTp1qnnxxRct36bffvutWbx4sUlNTTXGGHPkyBHTs2dP\n8+KLL5ri4mKTnZ1t+vfvbzZv3myMMSY1NdXMmTPHlJWVmbKyMrN161ZjjDHl5eXm//2//2feeOMN\nU1JSYo4dO2bGjh1rPvnkE2OMMV9//bV54IEHzP79+43L5TIff/yxGTZsmCktLTVlZWXmgQceMG+9\n9ZYpKyszGRkZpk+fPmb69OmW9v/s7TBu3Dhz7Ngxc/LkSTNy5EgzcuRI8+233xqXy2XS0tLMxIkT\nL9rH7Oxsk5GRYUpLS83JkyfNhAkTzMsvv3zO9j5bz549zeTJk83p06dNTk6OGTBggHnooYfMDz/8\nYMrKysxf//pXM2vWLGPMxT+fe/bsMdu3bzcul8vk5OSYMWPGmOXLl1doKysry+vbLysry/Tp08fz\n3sePHzeZmZnms88+Mw888IDJysoyTqfTzJkzx4wfP94YY0xeXp7p06ePycvL87zPgw8+aNauXWuM\n+en3tk+fPmbTpk2mrKzMLFu2zAwfPty4XC6v/BzOty+xwoW2zdm/c8YYU1JSYnr27GmOHDni6X//\n/v3Ntm3bjMvlMq+++qp5/PHHPc8/e99sjKnwfhfbJw4fPtyzz3vsscfM8OHDPZ/LSZMmmc8++8wY\nY8xrr71mJk2aZE6cOGGcTqeZPn26ef75540x/91PpKammsLCQuN0Oi3ZfhfC1ImFdu/erZycHA0b\nNkzh4eEKDAz0/DXdrVs3hYeHKygoSP3799fBgwd16tQpSVJgYKCys7N16tQphYWFqUmTJpJ++guz\nVatWuvXWWxUQEKAmTZro1ltv1bp1687bvt1u1+DBgxUUFKTGjRvruuuu8/zV4Au9evVSnTp1FBYW\npnbt2mnv3r364osvlJycrJtvvlkBAQG64447dN1112n9+vWe191888268cYbZbPZFBISosDAQOXm\n5iovL0/BwcFKSkqSJH3zzTeKiYnR3XffrcDAQNWtW1fdu3fXmjVrPO9Vv359denSRQEBAeratasc\nDoe++eYb5efna8OGDRoxYoTCw8MVGRmpfv36eV4bGRmp9u3bKyQkRGFhYerbt6+2bdtWoX9du3ZV\nQkKCgoKCdPvtt2vv3r0+2KrnN3DgQIWGhiouLk7NmzfXvn37JP30WcrLy1Nubq4CAwPVvHlzSdKe\nPXt07NgxDR48WMHBwapVq5Z69+7t6f/HH3+s++67T4mJiQoICNBdd90lm82mXbt2aefOnSopKVHf\nvn0VGBio1q1b+/zyyb///e9Vq1YtRUVFKTk5WbGxsWrTpo0CAgLUsWNH7d2796J9jIuLU+vWrRUU\nFKSoqCj16tXrnJ/x+fTp00eRkZFyOBxKSkpS48aN1ahRIwUGBqpDhw6ebX+xz2eDBg3UrFkzBQQE\nyOFw6Le//a22bt1q3Ub7PwEBAZJ+GnEsKSlRTEyM6tWrp1WrVqlXr16qW7euQkJCNHToUP3www/K\nzs5WzZo1ddNNN2n16tWSfvr85OXlqW3btp73bdOmjW666SYFBgaqb9++Kioq0q5du7zyczjfvsSX\n2+ZS3HLLLUpKSvLs1858Di7mYvvEli1basuWLXI6ncrOzlaPHj20detWuVwu7dixQy1atJAxRp98\n8omGDRum6OhohYSEaNCgQVq3bl2FaZIBAwYoPDxcISEhl7llfh2mTix07Ngx1alT55yh0PLyci1a\ntEjr1q3TyZMnPUPVp06dUlRUlJ544gm99dZbGjlypOLj4zVgwADdfPPNys3N1ddff60BAwZUeK/O\nnTuft/3IyMgKbYeEhMjpdHq/oxcQHR19TtvHjx8/Z3ohNjZWeXl5nsdnpifOGDt2rJYsWaKHH35Y\nUVFR6t27t+68807l5ORo3759FbaH2+32BJHzvZfD4fB88brdbj3wwAOedcYYud1uSVJJSYlee+01\nZWRkqKCgQJJUXFyssrIyzzatWbPmOf3zl5/XUlxcLEkaOnSolixZoqeeekoBAQH67W9/q759+yo3\nN1cnT57UwIEDPa9zu92eIf/c3FwtWLBAb775pmd9WVmZjh8/LpvNppiYGM/nVvppOxcWFlrdTY+f\nf7bO97O4WB9PnDih1157Tdu2bVNxcbGMMZd0o8Wft/XzWs5s+4t9Pg8dOqR58+Zpz549KikpUXl5\nuerXr3+5m+KyxcfHa9y4cfrwww81c+ZMNW/eXCkpKef8boaFhalGjRo6fvy44uPj1a1bNy1dulS9\ne/fWypUr1bFjxwr7l7N/1wICAlSrVi3P5+XX/hzOty+xwoW2zaWo7P7gYvvEFi1a6PPPP1fjxo3V\npEkTtWrVSrNnz9bu3bsVFRWl2NhY5efnq6SkROPHj6/wPjabTfn5+Z7HP98f+gpBw0K1a9fW0aNH\n5XK5FBj43029Zs0affXVV/rLX/6i2NhYFRUVVdgZNWjQQE888YTKy8u1du1a/e1vf9PixYtVp04d\ndezY8bxzfFeKWrVqnfPXSE5Ojlq1auV5/PMDpuLj4/XII4/IGKOtW7dq8uTJat68uerUqaOmTZvq\nueeeu2B7R48ePedxu3btVKdOHQUEBOiNN94475z4u+++q6ysLE2bNk0xMTHav3+/HnroIZkr7EK6\n0dHRnrOR9u3bp4kTJ6pRo0aqXbu2ateurddee+28r6tdu7buu+8+de3a9Zx1W7duVV5entxutyds\nHD169JLv5OgrF+vjm2++Kbfbrb///e+KiorS+vXr9dJLL3mt/Yt9Pl9++WXPsT3h4eH64IMPKozG\nWalDhw7q0KGDSkpKtHDhQs2ePVu1atVSbm6u5znFxcU6ffq05/ib3/zmN5ozZ4727t2rNWvW6Kmn\nnqrwnmf/rpWXl+v48eOqVauW7Ha7X38Ol+t82+aWW25RSUmJ5zknTpzwWnsX2ye2bNlSc+bM0aZN\nm9SyZUslJCTo6NGj2rBhg2eEPCoqSsHBwZoxY4ZiY2PPaSMnJ0fSuftWX2HqxEKNGjVSnTp1NG/e\nPBUVFcnlcmnr1q0qLi5WUFCQatSoodLSUi1atMjzmrKyMq1cuVIFBQUKCAhQRESEbDab7Ha7Onfu\nrI0bN2rDhg0qLy9XWVmZdu3apaysLD/28vJ06NBBGRkZ2rhxo8rLy7Vy5UodOnRIt9566wVfs2LF\nCuXn58tmsykiIkLST9NCZ0Z50tPTVVZWpvLych04cEC7d+/2vHb//v1avXq1p60jR47oN7/5jWrW\nrKnf/OY3mjt3rgoKCmSM0dGjRz3X/iguLlZwcLAiIiJUUFCgt956y9oNY5G1a9d6vgAiIiJkt9tl\nt9vVqFEjRUVFadmyZXI6nXK73Tp8+LBn6P53v/ud/vnPf2r//v0yxqi4uFgbNmxQUVGRmjZtquDg\nYC1fvlwul0ubN29WRkaGP7t5XhfrY3FxsUJDQxUREaHjx4/rvffeq/D66OhoHTlypNLtX+zzWVxc\nrPDwcIWFhenQoUP65JNPKt/Zy/Djjz9q8+bNKi0tVWBgoEJDQ2W329WpUyd98MEHOnTokEpLS7Vw\n4UI1atRI8fHxkn6ahuvcubNmzpypGjVqeKZ0z8jIyNB3330nl8uld955R2FhYZ5ppV/zc/ClC22b\nG264Qdu2bdORI0fkdDq1bNkyr7V5sX1i7dq1FRMTo/T0dLVs2VI2m01NmzbVxx9/rJYtW0r6aX94\n1113ad68eZ6RkPz8/ApT0v7EiIaFAgICNHHiRL322msaMWKE3G63WrZsqXHjxikjI0NDhw5VjRo1\nzrlt/erVqzV37lyVl5fL4XBo/PjxCg4OVu3atTV58mQtWLBAM2fOlCQlJiZWGP6v6q677jo9/vjj\nWrBggY4ePaprr71WEydOVFRU1AVfs3nzZr3++usqKSlRzZo1NWrUKMXFxUmS/vKXv2j+/Pl68803\n5XK5VLduXQ0ePNjz2rZt2+rbb7/VSy+9pNq1a2vChAmemxeNGzdOixYt0kMPPaSCggLVrl1bd911\nl9q0aaNevXrphRde0J/+9CfFxMSod+/e+vLLL63dOBbYs2eP5s2bp4KCAtWoUUM9e/b07JwmTpyo\nBQsWaNSoUXI6nYqNjVWfPn0kSbfeeqtKSko0Y8YM5eTkKCQkRElJSWrRooUCAwP15JNPavbs2frH\nP/6hFi1aqEuXLiorK/NnV89x5vfvQn0cMGCApk+frgEDBig+Pl6dO3fWu+++63n9mbMKFi1apB49\nelT4XF2KsLCwX/x8Dh06VHPmzNH777+vG264Qe3bt9emTZu8twEuoKysTIsWLVJWVpbnS3TMmDG6\n9tprdeLECU2ePNlz1snPh+K7du2q999/X/fff/8579u5c2e99957evbZZ1W3bl098cQTnpHcX/Nz\n8KULbZu6devqjjvu0MMPP6yIiAgNHjxYn3/+uVfavJR9YsuWLfXVV195zkw88/jMiIYkDRkyRG+/\n/bb+/Oc/6+TJk7rmmmvUsWPHX/wjzle4qRoA4JKcPn1aQ4YM0auvvuo5xgK4GKZOAAAXZYzRBx98\noOTkZEIGLgtTJwCAX1RWVqZBgwapZs2amjhxor/LwRWGqRMAAGAZpk4AAIBlCBoAAMAyBA0AAGAZ\nggYAALAMQQMAAFiGoAGg0qZPn6569eopICBAvXv3try93r17X/AmggCqJq6jAaBSdu/erUceeUSP\nP/64evbsyUWcAJwXQQNApezatUvGGA0fPlw33HCDv8sBUEUxdQLgsv3P//yPevbsKUlq0KCBbDab\nFixYoPz8fI0ZM0bx8fEKCQlRmzZtlJ6efs7rP/roI91yyy0KCwtTnTp1NHr0aBUWFlZ4zo4dO9Sp\nUyeFhoaqQYMGWrhwoU/6BsC7GNEAcNkmTpyopKQkPf7441q+fLni4+NVv3593XnnncrJydGzzz6r\n6667znPn04yMDM9dY//5z3/qj3/8o4YOHaopU6YoOztbf/7zn3XixAnP7bedTqe6d++uiIgIvfnm\nm5KkSZMm6dSpU2rUqJHf+g2gEgwAVMK7775rJJn9+/cbY4yZP3++CQwMNNu2bavwvFtuucX84Q9/\nMMYY43a7TUJCghkwYECF53z88cfGZrOZrVu3GmOMefnll43dbjc//PCD5zm7d+82drvddOrUybpO\nAfA6pk4AeEV6erpatmypxo0by+Vyef7deeed+uabbyRJP/zwgzIzM9WvX78Kz+nUqZPsdru+/fZb\nSdLXX3+tFi1aVBi9aNiwoVq1auWXvgGoPKZOAHjFsWPHtGnTJgUFBZ2zLiAgwPMcSbr33nvP+x5Z\nWVmSpOzsbDkcjnPWx8bGqri42FslA/ABggYAr4iJidGNN96oefPm/eJzJGn27Nm65ZZbzll/7bXX\nSpLi4+OVkZFxzvqcnBxFRUV5qWIAvkDQAOAV3bp107///W9de+21nsDwc02bNlXdunW1b98+Pfjg\ngxd8r7Zt2+qNN97Qnj171LBhQ0nSnj179N1336ljx46W1A/AGjZjjPF3EQCuPO+9957uvfde7d+/\nX4mJiSopKVH79u116tQpPfroo2rcuLHy8/O1adMmlZaWaurUqZKkt99+WwMHDtQDDzygHj16KCIi\nQpmZmfroo4/03HPPqXHjxiouLlbDhg1Vo0YNPfPMM5IqnnWyatUqP/YcwOVgRAOAV4SEhGjFihV6\n+umn9eyzzyo7O1u1a9dW69atNWbMGM/z/vCHPyg6OlrPPvusFi1aJElKTEzUXXfdpdjYWElSWFiY\n0tPTNXr0aA0ePFjXXXedJk6cqPfff1/5+fl+6R+AymFEAwAAWIbTWwEAgGUIGgAAwDIEDQAAYBmC\nBgAAsAxBAwAAWIagAQAALEPQAAAAliFoAAAAyxA0AACAZQgaAADAMgQNAABgmf8POkFhk0EPcNAA\nAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ggplot(chickens, aes(x=\"feed\",y=\"weight\")) + geom_boxplot() + theme_classic()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3. Null and alternative hyptheses for the difference in **chick weight** when fed **soybean feed** vs **sunflower feed**\n", + "\n", + "#### *Null hypothesis: mean(sunflower)-mean(soybean) = 0*\n", + "\n", + "#### *Alternative hypothesis: mean(sunflower)-mean(soybean) =/= 0*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 4. Test Null and alternative hypotheses using chi-2 likelihood test\n", + "\n", + "#### Create a data \"sun_soy_chickens\" containing only data of interest (chickens fed on soybean or sunflower feed). Add a column \"Factor\" to this new data frame, containing a 1 for chickens fed with sunflower feed, and a 0 for chickens fed with soybean feed." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tb/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: \n", + "A value is trying to be set on a copy of a slice from a DataFrame.\n", + "Try using .loc[row_indexer,col_indexer] = value instead\n", + "\n", + "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", + " This is separate from the ipykernel package so we can avoid doing imports until\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
weightfeedFactor
22243soybean0
23230soybean0
24248soybean0
25327soybean0
26329soybean0
27250soybean0
28193soybean0
29271soybean0
30316soybean0
31267soybean0
32199soybean0
33171soybean0
34158soybean0
35248soybean0
36423sunflower1
37340sunflower1
38392sunflower1
39339sunflower1
40341sunflower1
41226sunflower1
42320sunflower1
43295sunflower1
44334sunflower1
45322sunflower1
46297sunflower1
47318sunflower1
\n", + "
" + ], + "text/plain": [ + " weight feed Factor\n", + "22 243 soybean 0\n", + "23 230 soybean 0\n", + "24 248 soybean 0\n", + "25 327 soybean 0\n", + "26 329 soybean 0\n", + "27 250 soybean 0\n", + "28 193 soybean 0\n", + "29 271 soybean 0\n", + "30 316 soybean 0\n", + "31 267 soybean 0\n", + "32 199 soybean 0\n", + "33 171 soybean 0\n", + "34 158 soybean 0\n", + "35 248 soybean 0\n", + "36 423 sunflower 1\n", + "37 340 sunflower 1\n", + "38 392 sunflower 1\n", + "39 339 sunflower 1\n", + "40 341 sunflower 1\n", + "41 226 sunflower 1\n", + "42 320 sunflower 1\n", + "43 295 sunflower 1\n", + "44 334 sunflower 1\n", + "45 322 sunflower 1\n", + "46 297 sunflower 1\n", + "47 318 sunflower 1" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sun_soy_chickens=chickens[(chickens.feed == \"sunflower\") | (chickens.feed == \"soybean\")]\n", + "\n", + "sun_soy_chickens[\"Factor\"]=0\n", + "for i in range(0, len(sun_soy_chickens.feed)):\n", + " if sun_soy_chickens.feed.iloc[i] == \"sunflower\":\n", + " sun_soy_chickens.iloc[i,2]=1\n", + " else:\n", + " sun_soy_chickens.iloc[i,2]=0\n", + " \n", + "sun_soy_chickens" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Define null and alternative hypothesis functions" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def nllalt(p,obs):\n", + " B0=p[0]\n", + " sigma=p[1]\n", + " B1=p[2]\n", + " \n", + " expected=B0+B1*obs.Factor\n", + " nll=-1*norm(expected,sigma).logpdf(obs.weight).sum()\n", + " return nll\n", + " \n", + " \n", + "def nllnull(p,obs):\n", + " B0=p[0]\n", + " sigma=p[1]\n", + " \n", + " expected=B0\n", + " nll=-1*norm(expected,sigma).logpdf(obs.weight).sum()\n", + " return nll" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Make guess, use minimize function to find optimal coefficients for alternative and null hypothesis functions. T test to determine the significance of the factor (presence of sunflower seed) in determining the weight of chickens in this study" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Optimization terminated successfully.\n", + " Current function value: 138.469162\n", + " Iterations: 200\n", + " Function evaluations: 363\n", + "Optimization terminated successfully.\n", + " Current function value: 145.240592\n", + " Iterations: 138\n", + " Function evaluations: 261\n" + ] + }, + { + "data": { + "text/plain": [ + "0.00023317672869682671" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "guess=[1,1,1]\n", + "\n", + "fit_alt=minimize(nllalt, guess, method=\"Nelder-Mead\", options={'disp': True},args=sun_soy_chickens)\n", + "fit_null=minimize(nllnull, guess, method=\"Nelder-Mead\", options={'disp': True},args=sun_soy_chickens)\n", + "\n", + "D=2*(fit_null.fun-fit_alt.fun)\n", + "1-chi2.cdf(x=D, df=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Interpreting the result of the likelihood ratio test\n", + "\n", + "#### Since the P-value is \\<0.05, there is a significant difference in chicken weights between chickens fed on sunflower and soybean feed." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/12.1.py b/12.1.py new file mode 100644 index 0000000..de433ce --- /dev/null +++ b/12.1.py @@ -0,0 +1,108 @@ + +# coding: utf-8 + +# # General description of "chickwts.txt" dataset. +# +# +# +# ### Import packages and data + +# In[20]: + +import pandas as pd +import numpy as np +from scipy.stats import chi2 +from scipy.optimize import minimize +from scipy.stats import norm +from plotnine import * + +chickens=pd.read_table("chickwts.txt", delimiter=",") + + +# ### 1. Calculate mean weight of chickens based on feed source +# +# #### create an empty dictionary for "means". Fill by using a for loop to loop through each type of chicken feed, and add the mean weight of chickens for each feed to the dictionary entry for that feed. + +# In[21]: + +means=dict() + +for i in np.unique(chickens.feed): + means[i]=np.mean(chickens[chickens.feed == i][["weight"]])[0] + +means_df=pd.DataFrame(means, index =[0]) +means_df + + +# ### 2. Generate a plot summarizing chicken weight data +# +# #### use ggplot2 to create a boxplot of "chickwts.txt" data, plotting chicken **weight** on the y axis, grouped by chicken **feed** + +# In[22]: + +ggplot(chickens, aes(x="feed",y="weight")) + geom_boxplot() + theme_classic() + + +# ### 3. Null and alternative hyptheses for the difference in **chick weight** when fed **soybean feed** vs **sunflower feed** +# +# #### *Null hypothesis: mean(sunflower)-mean(soybean) = 0* +# +# #### *Alternative hypothesis: mean(sunflower)-mean(soybean) =/= 0* + +# ### 4. Test Null and alternative hypotheses using chi-2 likelihood test +# +# #### Create a data "sun_soy_chickens" containing only data of interest (chickens fed on soybean or sunflower feed). Add a column "Factor" to this new data frame, containing a 1 for chickens fed with sunflower feed, and a 0 for chickens fed with soybean feed. + +# In[23]: + +sun_soy_chickens=chickens[(chickens.feed == "sunflower") | (chickens.feed == "soybean")] + +sun_soy_chickens["Factor"]=0 +for i in range(0, len(sun_soy_chickens.feed)): + if sun_soy_chickens.feed.iloc[i] == "sunflower": + sun_soy_chickens.iloc[i,2]=1 + else: + sun_soy_chickens.iloc[i,2]=0 + +sun_soy_chickens + + +# #### Define null and alternative hypothesis functions + +# In[24]: + +def nllalt(p,obs): + B0=p[0] + sigma=p[1] + B1=p[2] + + expected=B0+B1*obs.Factor + nll=-1*norm(expected,sigma).logpdf(obs.weight).sum() + return nll + + +def nllnull(p,obs): + B0=p[0] + sigma=p[1] + + expected=B0 + nll=-1*norm(expected,sigma).logpdf(obs.weight).sum() + return nll + + +# #### Make guess, use minimize function to find optimal coefficients for alternative and null hypothesis functions. T test to determine the significance of the factor (presence of sunflower seed) in determining the weight of chickens in this study + +# In[25]: + +guess=[1,1,1] + +fit_alt=minimize(nllalt, guess, method="Nelder-Mead", options={'disp': True},args=sun_soy_chickens) +fit_null=minimize(nllnull, guess, method="Nelder-Mead", options={'disp': True},args=sun_soy_chickens) + +D=2*(fit_null.fun-fit_alt.fun) +1-chi2.cdf(x=D, df=1) + + +# ### Interpreting the result of the likelihood ratio test +# +# #### Since the P-value is \<0.05, there is a significant difference in chicken weights between chickens fed on sunflower and soybean feed.