<jupyter_start><jupyter_text>Statistics There are many specialized packages for dealing with data analysis and statistical programming. One very important code that you will see in MATH1024, Introduction to Probability and Statistics, is [R](http://www.r-project.org/). A Python package for performing similar analysis of large data sets is [pandas](http://pandas.pydata.org/). However, simple statistical tasks on simple data sets can be tackled using `numpy` and `scipy`.  Getting data in A data file containing the monthly rainfall for Southampton, taken from the [Met Office data](http://www.metoffice.gov.uk/pub/data/weather/uk/climate/stationdata/southamptondata.txt) can be [downloaded from this link](https://github.com/IanHawke/maths-with-python/blob/master/southampton_precip.txt). We will save that file locally, and then look at the data.The first few lines of the file are:<jupyter_code>!head southampton_precip.txt<jupyter_output>#Year Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
1855 85.6 54.3 61.3 10.1 60.0 43.9 101.0 47.9 88.4 187.5 28.2 55.4
1856 93.5 50.6 36.3 127.3 55.7 40.3 16.5 64.7 67.6 74.5 38.7 87.1
1857 72.3 10.6 54.4 60.7 19.0 38.2 43.7 66.3 93.6 191.4 57.1 25.0
1858 27.0 33.1 22.9 94.1 65.7 14.1 69.6 55.5 75.2 66.2 50.1 116.6
1859 59.6 78.3 49.7 92.4 36.8 45.7 66.6 58.3 135.3 119.8 125.1 127.1
1860 129.2 29.3 59.3 47.6 88.7 205.0 84.7 115.0 99.2 53.2 80.2 127.7
1861 20.7 60.2 76.4 10.2 41.3 100.8 103.5 22.2 78.0 27.7 164.3 53.2
1862 104.0 20.1 124.2 57.5 123.9 53.8 52.8 36.3 29.7 171.8 22.4 72.7
1863 129.4 32.4 38.7 20.5 55.2 94.6 26.4 63.9 98.7 115.3 60.7 64.4<jupyter_text>We can use `numpy` to load this data into a variable, where we can manipulate it. This is not ideal: it will lose the information in the header, and that the first column corresponds to years. However, it is simple to use.<jupyter_code>import numpy
data = numpy.loadtxt('southampton_precip.txt')
data<jupyter_output><empty_output><jupyter_text>We see that the first column - the year - has been converted to a floating point number, which is not helpful. However, we can now split the data using standard `numpy` operations:<jupyter_code>years = data[:, 0]
rainfall = data[:, 1:]<jupyter_output><empty_output><jupyter_text>We can now plot, for example, the rainfall in January for all years:<jupyter_code>%matplotlib inline
from matplotlib import rcParams
rcParams['figure.figsize']=(12,9)
from matplotlib import pyplot
pyplot.plot(years, rainfall[:,0])
pyplot.xlabel('Year')
pyplot.ylabel('Rainfall in January');<jupyter_output><empty_output><jupyter_text>Basic statistical functions `numpy` contains a number of basic statistical functions, such as `min`, `max` and `mean`. These will act on entire arrays to give the "all time" minimum, maximum, and average rainfall:<jupyter_code>print("Minimum rainfall: {}".format(rainfall.min()))
print("Maximum rainfall: {}".format(rainfall.max()))
print("Mean rainfall: {}".format(rainfall.mean()))<jupyter_output>Minimum rainfall: 0.0
Maximum rainfall: 280.7
Mean rainfall: 67.03591954022988<jupyter_text>Of more interest would be either1. the mean (min/max) rainfall in a given month for all years, or2. the mean (min/max) rainfall in a given year for all months.So the mean rainfall in the first year, 1855, would be<jupyter_code>print ("Mean rainfall in 1855: {}".format(rainfall[0, :].mean()))<jupyter_output>Mean rainfall in 1855: 68.63333333333334<jupyter_text>Whilst the mean rainfall in January, averaging over all years, would be<jupyter_code>print ("Mean rainfall in January: {}".format(rainfall[:, 0].mean()))<jupyter_output>Mean rainfall in January: 81.86482758620689<jupyter_text>If we wanted to plot the mean rainfall per year, across all years, this would be tedious - there are 145 years of data in the file. Even computing the mean rainfall in each month, across all years, would be bad with 12 months. We could write a loop. However, `numpy` allows us to apply a function along an axis of the array, which does this is one operation:<jupyter_code>mean_rainfall_in_month = rainfall.mean(axis=0)
mean_rainfall_per_year = rainfall.mean(axis=1)<jupyter_output><empty_output><jupyter_text>The `axis` argument gives the direction we want to *keep* - that we do not apply the operation to. For this data set, each row contains a year and each column a month. To find the mean in a given month we want to keep the row information (`axis` 0) and take the mean over the column. To find the mean in a given year we want to keep the column information (`axis` 1) and take the mean over the row.We can now plot how the mean varies with each year.<jupyter_code>pyplot.plot(years, mean_rainfall_per_year)
pyplot.xlabel('Year')
pyplot.ylabel('Mean rainfall');<jupyter_output><empty_output><jupyter_text>We can also compute the standard deviation:<jupyter_code>std_rainfall_per_year = rainfall.std(axis=1)<jupyter_output><empty_output><jupyter_text>We can then add confidence intervals to the plot:<jupyter_code>pyplot.errorbar(years, mean_rainfall_per_year, yerr = std_rainfall_per_year)
pyplot.xlabel('Year')
pyplot.ylabel('Mean rainfall');<jupyter_output><empty_output><jupyter_text>This isn't particularly pretty or clear: a nicer example would use better packages, but a quick fix uses an alternative `matplotlib` approach:<jupyter_code>pyplot.plot(years, mean_rainfall_per_year)
pyplot.fill_between(years, mean_rainfall_per_year - std_rainfall_per_year, 
                    mean_rainfall_per_year + std_rainfall_per_year,
                    alpha=0.25, color=None)
pyplot.xlabel('Year')
pyplot.ylabel('Mean rainfall');<jupyter_output><empty_output><jupyter_text>Categorical data Looking at the means by month, it would be better to give them names rather than numbers. We will also summarize the available information using a boxplot:<jupyter_code>months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
pyplot.boxplot(rainfall, labels=months)
pyplot.xlabel('Month')
pyplot.ylabel('Mean rainfall');<jupyter_output><empty_output><jupyter_text>Much better ways of working with categorical data are available through more specialized packages.  Regression We can go beyond the basic statistical functions in `numpy` and look at other standard tasks. For example, we can look for simple trends in our data with a *linear regression*. There is a function to compute the linear regression in `scipy` we can use. We will use this to see if there is a trend in the mean yearly rainfall:<jupyter_code>from scipy import stats

slope, intercept, r_value, p_value, std_err = stats.linregress(years, mean_rainfall_per_year)

pyplot.plot(years, mean_rainfall_per_year, 'b-', label='Data')
pyplot.plot(years, intercept + slope*years, 'k-', label='Linear Regression')
pyplot.xlabel('Year')
pyplot.ylabel('Mean rainfall')
pyplot.legend();
print("The change in rainfall (the slope) is {}.".format(slope))
print("However, the error estimate is {}.".format(std_err))
print("The correlation coefficient between rainfall and year"
      " is {}.".format(r_value))
print("The probability that the slope is zero is {}.".format(p_value))<jupyter_output>The change in rainfall (the slope) is -0.028739338949246847.
However, the error estimate is 0.021587122926201515.
The correlation coefficient between rainfall and year is -0.11064686384415015.
The probability that the slope is zero is 0.18520267346715713.<jupyter_text>It looks like there's a good chance that the slight decrease in mean rainfall with time is a real effect.  Random numbers Random processes and random variables may be at the heart of probability and statistics, but computers cannot generate anything "truly" random. Instead they can generate *pseudo-random* numbers using random number generators (RNGs). Constructing a random number generator is a *hard problem* and wherever possible you should use a well-tested RNG rather than attempting to write your own.Python has many ways of generating random numbers. Perhaps the most useful are given by the [`numpy.random`](http://docs.scipy.org/doc/numpy/reference/routines.random.html) module, which can generate a `numpy` array filled with random numbers from various distributions. For example:<jupyter_code>from numpy import random

uniform = random.rand(10000)
normal = random.randn(10000)

fig = pyplot.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.hist(uniform, 20)
ax1.set_title('Uniform data')
ax2.hist(normal, 20)
ax2.set_title('Normal data')
fig.tight_layout()
fig.show();<jupyter_output>/Users/ih3/anaconda/lib/python3.4/site-packages/matplotlib/figure.py:397: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "<jupyter_text>More distributions Whilst the standard distributions are given by the convenience functions above, the [full documentation of `numpy.random`](http://docs.scipy.org/doc/numpy/reference/routines.random.html) shows many other distributions available. For example, we can draw $10,000$ samples from the [Beta distribution](http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.beta.htmlnumpy.random.beta) using the parameters $\alpha = 1/2 = \beta$ as<jupyter_code>beta_samples = random.beta(0.5, 0.5, 10000)

pyplot.hist(beta_samples, 20)
pyplot.title('Beta data')
pyplot.show();<jupyter_output><empty_output><jupyter_text>We can do this $5,000$ times and compute the mean of each set of samples:<jupyter_code>n_trials = 5000
beta_means = numpy.zeros((n_trials,))

for trial in range(n_trials):
    beta_samples = random.beta(0.5, 0.5, 10000)
    beta_means[trial] = numpy.mean(beta_samples)
    
pyplot.hist(beta_means, 20)
pyplot.title('Mean of Beta trials')
pyplot.show();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1021763563911311, 'fraction_numerical': 0.06314498824971901, 'mean_word_length': 5.258312020460358, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 63, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '74453', 'license': 'MIT', 'n_tokens_mistral': 3642, 'n_tokens_neox': 3122, 'n_words': 1118, 'path': '08-statistics.ipynb', 'repo_name': 'LaGuer/maths-with-python', 'size': '9787'}