#!/usr/bin/env python

# ----------------------------------------------------------------------------------- #
#
#  Python macro for exercise on planning experiment.
#
#  You are planning the "Lady Tasting Tea" experiment, serving more cups of tea (or
#  whatever you're tasting/testing). However, you are wondering, if you experiment is
#  good enough for the task defined, and would like to improve on it (also beyond
#  simply tasting tea).
#
#  In this exercise, you are welcome to drop the idea of "tasting/testing", and make
#  estimates for YOUR favorite dream experiment. I simply want you to try to ask some
#  questions about the performance of your experiment before doing it, and then testing
#  which requirements this puts on your experimental plan/design. But if you don't have
#  any ideas for this, you can just follow the exercise below.
#
#  Author: Troels C. Petersen (NBI)
#  Email:  petersen@nbi.dk
#  Date:   2nd of January 2018
#
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
# Imports
# ----------------------------------------------------------------------------------- #

# First, import the modules you want to use:
from __future__ import print_function, division        #
import numpy as np                                     # Matlab like syntax for linear algebra and functions
import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab
from scipy.special import gamma, binom
plt.close('all')         # to close all open figures

r = np.random             # Random generator
r.seed(42)                # Set a random seed (but a fixed one)
SavePlots = True          # For now, don't save plots (once you trust your code, switch on)


Nexperiments = 100
Ntests = 8
Ncorr  = 4


# 10-cup lady-tasting-tea experiment:
# -----------------------------------
# With two options:
#   Binomial: N=10, p=0.5
#


Npos = int(gamma(Ntests+1) / gamma(Ncorr+1) / gamma(Ntests-Ncorr+1))
for i in range (Ntests) :
    Npos_i = binom(Ntests,i)
    print(i, Npos_i, Npos, Npos_i/Npos)

#
# If more options are involved:
#   Multinomial: N=10, p1=1/3 p2=1/3 p3=1/3



# ------------------------------------------------------------------ #
# Loop over number of experiments and make random choices:
# ------------------------------------------------------------------ #

Ncorrect_all = []

for iexp in range (Nexperiments) :
    
    # Make Ntests tests of a random number:
    Ncorrect = 0
    for itest in range (Ntests) :
        if (r.uniform() > 0.5) : 
            Ncorrect += 1
            
    Ncorrect_all.append(Ncorrect)
#    Hist_Ncorrect.Fill(Ncorrect)
    
Ncorrect_all = np.array(Ncorrect_all)



# ------------------------------------------------------------------ #
# Show the distribution of results:
# ------------------------------------------------------------------ #

fig, ax = plt.subplots()
ax.hist(Ncorrect_all, Ntests+1, (-0.5, Ntests+0.5), histtype='step')
ax.set(xlabel='Number of correct answers', ylabel='Frequency', title='Histogram of Ncorrect')

fig.tight_layout()
plt.show(block=False)

if SavePlots:
    fig.savefig('TestResult.pdf', dpi=600)

try:
    __IPYTHON__
except:
    raw_input('Press Enter to exit')



# ----------------------------------------------------------------------------------- #
# 
#  Before starting looking at the program, sit down and calculate (either in hand, or
#  as a section in the beginning of this program) the probabilities of the outcome of
#  a 10-cup lady-tasting-tea experiment. Such a calculation is always healthy to get
#  a feel for the sensitivity of the experiment.
#
#  Questions:
#  ----------
#   1) What is the chance that a person choosing at random gets all 10 tastings correct?
#      And given that this actually happens once in the above program, would you claim,
#      that I chose the "random" seed deliberately and not at random?
#
#   2) In the above example, each test is random, and there is no knowledge of how many
#      of each sample is being tested. How does the test change, when it is known, that
#      half of the samples are of one type and the other half of the opposite type?
# 
#   3) Imagine that you are now conducting a beer tasting! You friend made the (drunk?)
#      claim, that he could always recognize a "Groen Tuborg". You decide to design an
#      experiment that can test this hypothesis. How would you do this?
#       - How sure do you want to be on the result? (i.e. p-value of all correct)
#       - How many types of beer would you use?
#       - What if you want to allow for one (or more) wrong identifications?
#       - How many times should your friend taste a beer sample?
#       - Will you tell him ahead of time, how many of each sample there are?
# 
#   4) Imagine a prominent professor at NBI making a similar claim about sushi? You
#      can only ask him to taste 12 times (that is a full lunch, and you can not
#      afford more). If you get three samples, four from each restaurant, and you
#      allow for up to one misidentification, how sure can you be, that his claim is
#      reasonable? I.e. what is the p-value for a random guess to get 11 or 12 correct?
#      Calculate this both analytically and using a Monte Carlo simulation.
# 
# ----------------------------------------------------------------------------------- #