!head southampton_precip.txt

"""
Explanation: 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. A Python package for performing similar analysis of large data sets is pandas. 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 can be downloaded from this link. We will save that file locally, and then look at the data.
The first few lines of the file are:
End of explanation
"""


import numpy

data = numpy.loadtxt('southampton_precip.txt')

data

"""
Explanation: 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.
End of explanation
"""


years = data[:, 0]
rainfall = data[:, 1:]

"""
Explanation: 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:
End of explanation
"""


%matplotlib inline
from matplotlib import pyplot

pyplot.plot(years, rainfall[:,0])
pyplot.xlabel('Year')
pyplot.ylabel('Rainfall in January');

"""
Explanation: We can now plot, for example, the rainfall in January for all years:
End of explanation
"""


print("Minimum rainfall: {}".format(rainfall.min()))
print("Maximum rainfall: {}".format(rainfall.max()))
print("Mean rainfall: {}".format(rainfall.mean()))

"""
Explanation: 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:
End of explanation
"""


print ("Mean rainfall in 1855: {}".format(rainfall[0, :].mean()))

"""
Explanation: Of more interest would be either

the mean (min/max) rainfall in a given month for all years, or
the mean (min/max) rainfall in a given year for all months.

So the mean rainfall in the first year, 1855, would be
End of explanation
"""


print ("Mean rainfall in January: {}".format(rainfall[:, 0].mean()))

"""
Explanation: Whilst the mean rainfall in January, averaging over all years, would be
End of explanation
"""


mean_rainfall_in_month = rainfall.mean(axis=0)
mean_rainfall_per_year = rainfall.mean(axis=1)

"""
Explanation: 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:
End of explanation
"""


pyplot.plot(years, mean_rainfall_per_year)
pyplot.xlabel('Year')
pyplot.ylabel('Mean rainfall');

"""
Explanation: 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.
End of explanation
"""


std_rainfall_per_year = rainfall.std(axis=1)

"""
Explanation: We can also compute the standard deviation:
End of explanation
"""


pyplot.errorbar(years, mean_rainfall_per_year, yerr = std_rainfall_per_year)
pyplot.xlabel('Year')
pyplot.ylabel('Mean rainfall');

"""
Explanation: We can then add confidence intervals to the plot:
End of explanation
"""


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');

"""
Explanation: This isn't particularly pretty or clear: a nicer example would use better packages, but a quick fix uses an alternative matplotlib approach:
End of explanation
"""


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');

"""
Explanation: 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:
End of explanation
"""


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 probability that the rainfall is correlated"
      " with the year is {}.".format(r_value))
print("The probability that the slope is zero is {}.".format(p_value))

"""
Explanation: 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:
End of explanation
"""


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();

"""
Explanation: It looks like there's no significant change in mean rainfall with time.
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 module, which can generate a numpy array filled with random numbers from various distributions. For example:
End of explanation
"""


beta_samples = random.beta(0.5, 0.5, 10000)

pyplot.hist(beta_samples, 20)
pyplot.title('Beta data')
pyplot.show();

"""
Explanation: More distributions
Whilst the standard distributions are given by the convenience functions above, the full documentation of numpy.random shows many other distributions available. For example, we can draw $10,000$ samples from the Beta distribution using the parameters $\alpha = 1/2 = \beta$ as
End of explanation
"""


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();

"""
Explanation: We can do this $5,000$ times and compute the mean of each set of samples:
End of explanation
"""


