from __future__ import print_function, division

import numpy
import scipy.stats

import matplotlib.pyplot as pyplot

from ipywidgets import interact, interactive, fixed
import ipywidgets as widgets

# seed the random number generator so we all get the same results
numpy.random.seed(18)

# some nicer colors from http://colorbrewer2.org/
COLOR1 = '#7fc97f'
COLOR2 = '#beaed4'
COLOR3 = '#fdc086'
COLOR4 = '#ffff99'
COLOR5 = '#386cb0'

%matplotlib inline

"""
Explanation: Random Sampling
Copyright 2016 Allen Downey
License: Creative Commons Attribution 4.0 International
End of explanation
"""


weight = scipy.stats.lognorm(0.23, 0, 70.8)
weight.mean(), weight.std()

"""
Explanation: Part One
Suppose we want to estimate the average weight of men and women in the U.S.
And we want to quantify the uncertainty of the estimate.
One approach is to simulate many experiments and see how much the results vary from one experiment to the next.
I'll start with the unrealistic assumption that we know the actual distribution of weights in the population.  Then I'll show how to solve the problem without that assumption.
Based on data from the BRFSS, I found that the distribution of weight in kg for women in the U.S. is well modeled by a lognormal distribution with the following parameters:
End of explanation
"""


xs = numpy.linspace(20, 160, 100)
ys = weight.pdf(xs)
pyplot.plot(xs, ys, linewidth=4, color=COLOR1)
pyplot.xlabel('weight (kg)')
pyplot.ylabel('PDF')
None

"""
Explanation: Here's what that distribution looks like:
End of explanation
"""


def make_sample(n=100):
    sample = weight.rvs(n)
    return sample

"""
Explanation: make_sample draws a random sample from this distribution.  The result is a NumPy array.
End of explanation
"""


sample = make_sample(n=100)
sample.mean(), sample.std()

"""
Explanation: Here's an example with n=100.  The mean and std of the sample are close to the mean and std of the population, but not exact.
End of explanation
"""


def sample_stat(sample):
    return sample.mean()

"""
Explanation: We want to estimate the average weight in the population, so the "sample statistic" we'll use is the mean:
End of explanation
"""


def compute_sampling_distribution(n=100, iters=1000):
    stats = [sample_stat(make_sample(n)) for i in range(iters)]
    return numpy.array(stats)

"""
Explanation: One iteration of "the experiment" is to collect a sample of 100 women and compute their average weight.
We can simulate running this experiment many times, and collect a list of sample statistics.  The result is a NumPy array.
End of explanation
"""


sample_means = compute_sampling_distribution(n=100, iters=1000)

"""
Explanation: The next line runs the simulation 1000 times and puts the results in
sample_means:
End of explanation
"""


pyplot.hist(sample_means, color=COLOR5)
pyplot.xlabel('sample mean (n=100)')
pyplot.ylabel('count')
None

"""
Explanation: Let's look at the distribution of the sample means.  This distribution shows how much the results vary from one experiment to the next.
Remember that this distribution is not the same as the distribution of weight in the population.  This is the distribution of results across repeated imaginary experiments.
End of explanation
"""


sample_means.mean()

"""
Explanation: The mean of the sample means is close to the actual population mean, which is nice, but not actually the important part.
End of explanation
"""


std_err = sample_means.std()
std_err

"""
Explanation: The standard deviation of the sample means quantifies the variability from one experiment to the next, and reflects the precision of the estimate.
This quantity is called the "standard error".
End of explanation
"""


conf_int = numpy.percentile(sample_means, [5, 95])
conf_int

"""
Explanation: We can also use the distribution of sample means to compute a "90% confidence interval", which contains 90% of the experimental results:
End of explanation
"""


def plot_sampling_distribution(n, xlim=None):
    """Plot the sampling distribution.
    
    n: sample size
    xlim: [xmin, xmax] range for the x axis 
    """
    sample_stats = compute_sampling_distribution(n, iters=1000)
    se = numpy.std(sample_stats)
    ci = numpy.percentile(sample_stats, [5, 95])
    
    pyplot.hist(sample_stats, color=COLOR2)
    pyplot.xlabel('sample statistic')
    pyplot.xlim(xlim)
    text(0.03, 0.95, 'CI [%0.2f %0.2f]' % tuple(ci))
    text(0.03, 0.85, 'SE %0.2f' % se)
    pyplot.show()
    
def text(x, y, s):
    """Plot a string at a given location in axis coordinates.
    
    x: coordinate
    y: coordinate
    s: string
    """
    ax = pyplot.gca()
    pyplot.text(x, y, s,
                horizontalalignment='left',
                verticalalignment='top',
                transform=ax.transAxes)

"""
Explanation: Now we'd like to see what happens as we vary the sample size, n.  The following function takes n, runs 1000 simulated experiments, and summarizes the results.
End of explanation
"""


plot_sampling_distribution(100)

"""
Explanation: Here's a test run with n=100:
End of explanation
"""


def sample_stat(sample):
    return sample.mean()

slider = widgets.IntSlider(min=10, max=1000, value=100)
interact(plot_sampling_distribution, n=slider, xlim=fixed([55, 95]))
None

"""
Explanation: Now we can use interact to run plot_sampling_distribution with different values of n.  Note: xlim sets the limits of the x-axis so the figure doesn't get rescaled as we vary n.
End of explanation
"""


def sample_stat(sample):
    # TODO: replace the following line with another sample statistic
    return sample.mean()

slider = widgets.IntSlider(min=10, max=1000, value=100)
interact(plot_sampling_distribution, n=slider, xlim=fixed([0, 100]))
None

"""
Explanation: Other sample statistics
This framework works with any other quantity we want to estimate.  By changing sample_stat, you can compute the SE and CI for any sample statistic.
Exercise 1: Fill in sample_stat below with any of these statistics:

Standard deviation of the sample.
Coefficient of variation, which is the sample standard deviation divided by the sample standard mean.
Min or Max
Median (which is the 50th percentile)
10th or 90th percentile.
Interquartile range (IQR), which is the difference between the 75th and 25th percentiles.

NumPy array methods you might find useful include std, min, max, and percentile.
Depending on the results, you might want to adjust xlim.
End of explanation
"""


class Resampler(object):
    """Represents a framework for computing sampling distributions."""
    
    def __init__(self, sample, xlim=None):
        """Stores the actual sample."""
        self.sample = sample
        self.n = len(sample)
        self.xlim = xlim
        
    def resample(self):
        """Generates a new sample by choosing from the original
        sample with replacement.
        """
        new_sample = numpy.random.choice(self.sample, self.n, replace=True)
        return new_sample
    
    def sample_stat(self, sample):
        """Computes a sample statistic using the original sample or a
        simulated sample.
        """
        return sample.mean()
    
    def compute_sampling_distribution(self, iters=1000):
        """Simulates many experiments and collects the resulting sample
        statistics.
        """
        stats = [self.sample_stat(self.resample()) for i in range(iters)]
        return numpy.array(stats)
    
    def plot_sampling_distribution(self):
        """Plots the sampling distribution."""
        sample_stats = self.compute_sampling_distribution()
        se = sample_stats.std()
        ci = numpy.percentile(sample_stats, [5, 95])
    
        pyplot.hist(sample_stats, color=COLOR2)
        pyplot.xlabel('sample statistic')
        pyplot.xlim(self.xlim)
        text(0.03, 0.95, 'CI [%0.2f %0.2f]' % tuple(ci))
        text(0.03, 0.85, 'SE %0.2f' % se)
        pyplot.show()

"""
Explanation: STOP HERE
We will regroup and discuss before going on.
Part Two
So far we have shown that if we know the actual distribution of the population, we can compute the sampling distribution for any sample statistic, and from that we can compute SE and CI.
But in real life we don't know the actual distribution of the population.  If we did, we wouldn't be doing statistical inference in the first place!
In real life, we use the sample to build a model of the population distribution, then use the model to generate the sampling distribution.  A simple and popular way to do that is "resampling," which means we use the sample itself as a model of the population distribution and draw samples from it.
Before we go on, I want to collect some of the code from Part One and organize it as a class.  This class represents a framework for computing sampling distributions.
End of explanation
"""


def interact_func(n, xlim):
    sample = weight.rvs(n)
    resampler = Resampler(sample, xlim=xlim)
    resampler.plot_sampling_distribution()

"""
Explanation: The following function instantiates a Resampler and runs it.
End of explanation
"""


interact_func(n=100, xlim=[50, 100])

"""
Explanation: Here's a test run with n=100
End of explanation
"""


slider = widgets.IntSlider(min=10, max=1000, value=100)
interact(interact_func, n=slider, xlim=fixed([50, 100]))
None

"""
Explanation: Now we can use interact_func in an interaction:
End of explanation
"""


# Solution goes here

"""
Explanation: Exercise 2: write a new class called StdResampler that inherits from Resampler and overrides sample_stat so it computes the standard deviation of the resampled data.
End of explanation
"""


def interact_func2(n, xlim):
    sample = weight.rvs(n)
    resampler = StdResampler(sample, xlim=xlim)
    resampler.plot_sampling_distribution()
    
interact_func2(n=100, xlim=[0, 100])

"""
Explanation: Test your code using the cell below:
End of explanation
"""


slider = widgets.IntSlider(min=10, max=1000, value=100)
interact(interact_func2, n=slider, xlim=fixed([0, 100]))
None

"""
Explanation: When your StdResampler is working, you should be able to interact with it:
End of explanation
"""


female_weight = scipy.stats.lognorm(0.23, 0, 70.8)
female_weight.mean(), female_weight.std()

"""
Explanation: STOP HERE
We will regroup and discuss before going on.
Part Three
We can extend this framework to compute SE and CI for a difference in means.
For example, men are heavier than women on average.  Here's the women's distribution again (from BRFSS data):
End of explanation
"""


male_weight = scipy.stats.lognorm(0.20, 0, 87.3)
male_weight.mean(), male_weight.std()

"""
Explanation: And here's the men's distribution:
End of explanation
"""


female_sample = female_weight.rvs(100)
male_sample = male_weight.rvs(100)

"""
Explanation: I'll simulate a sample of 100 men and 100 women:
End of explanation
"""


male_sample.mean() - female_sample.mean()

"""
Explanation: The difference in means should be about 17 kg, but will vary from one random sample to the next:
End of explanation
"""


def CohenEffectSize(group1, group2):
    """Compute Cohen's d.

    group1: Series or NumPy array
    group2: Series or NumPy array

    returns: float
    """
    diff = group1.mean() - group2.mean()

    n1, n2 = len(group1), len(group2)
    var1 = group1.var()
    var2 = group2.var()

    pooled_var = (n1 * var1 + n2 * var2) / (n1 + n2)
    d = diff / numpy.sqrt(pooled_var)
    return d

"""
Explanation: Here's the function that computes Cohen's effect size again:
End of explanation
"""


CohenEffectSize(male_sample, female_sample)

"""
Explanation: The difference in weight between men and women is about 1 standard deviation:
End of explanation
"""


class CohenResampler(Resampler):
    def __init__(self, group1, group2, xlim=None):
        self.group1 = group1
        self.group2 = group2
        self.xlim = xlim
        
    def resample(self):
        n, m = len(self.group1), len(self.group2)
        group1 = numpy.random.choice(self.group1, n, replace=True)
        group2 = numpy.random.choice(self.group2, m, replace=True)
        return group1, group2
    
    def sample_stat(self, groups):
        group1, group2 = groups
        return CohenEffectSize(group1, group2)

"""
Explanation: Now we can write a version of the Resampler that computes the sampling distribution of $d$.
End of explanation
"""


resampler = CohenResampler(male_sample, female_sample)
resampler.plot_sampling_distribution()

"""
Explanation: Now we can instantiate a CohenResampler and plot the sampling distribution.
End of explanation
"""


