<jupyter_start><jupyter_text>Descriptive statistics<jupyter_code>import numpy as np
import seaborn as sns
import scipy.stats as st
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import pandas as pd
import statsmodels.api as sm
import statistics
import os 
from scipy.stats import norm<jupyter_output><empty_output><jupyter_text>Probability data, binomial distribution We already got to know data that follow a binomial distribution, but we actually had not looked at the distribution. We will do this now. 10% of the 100 cells we count have deformed nuclei. To illustrate the distribution we will count repeatedly....<jupyter_code>n = 100 # number of trials
p = 0.1 # probability of each trial

s = np.random.binomial(n, p, 1000) #simulation repeating the experiment 1000 times
print(s)<jupyter_output>[ 8  6  9  9 12  8 14 10 18 11 10 10 13 11 12 10 11 12 10  5  6 12 13 12
  8  9  9  6 14 12  5 11  6 10 10  9 16 11  9 15 13 12 10  9 10 10 10  9
 11  9 10 10  9  8  9 13 10 10 12  8 11  9  6 13  6 12  7  9  7 13 10  7
 10 12 21 10 11  9 11  8  8 11  8 11 11 19  8 13 11 11 10 16 10  9  6 11
  9  7 13  5 12  7  6 13 10 15  8  6 10  8 12  8 10 10 11 12 11 12  8  8
 11  8 14 11 10 10 10  7  5 13 14  9  6 13 10  7 18 12 14 12  9 13 10 13
  8 10  9 10  9  9  8  7  6 11  8 11 10 10  7 10 21 13 16 14 10 13  9 15
 13  7 13  6  5  8 12 11  8 11  8 11  6 12  7 11 14 13 15  6 13 11 11 11
 15 13 14  5 11  9  8  6 14  6  8  8  8 14 12  6 14 10 13 11 16  7  8 14
  9  7  4 11 12 14 11 12  7 13  5  8  9  9  4 13 12  6 13  7 15 11 10  9
 15 14  9 17 10 11  9  8 15 12 12 10 12 15  8 14  9 11  9 15  9 15 12  9
  5 11  9  7  9 14  8  9  7 11  8  6 11 13  6 16 14 12  7 15  8  9 10 10
  9 14 11 10 16 10  9 11 12 14 18 10  5  9  9 11 12 10 10 12  9 10 11 10
 12 10 10  9  9  5  5  7 13  7 12 10  6 10  8  9  8[...]<jupyter_text>As you can see, the result of the distribution is in absolute counts, not proportions - they can easyly be converted by deviding with n, but they dont have to...<jupyter_code>props = s/n
print(props)<jupyter_output>[0.08 0.06 0.09 0.09 0.12 0.08 0.14 0.1  0.18 0.11 0.1  0.1  0.13 0.11
 0.12 0.1  0.11 0.12 0.1  0.05 0.06 0.12 0.13 0.12 0.08 0.09 0.09 0.06
 0.14 0.12 0.05 0.11 0.06 0.1  0.1  0.09 0.16 0.11 0.09 0.15 0.13 0.12
 0.1  0.09 0.1  0.1  0.1  0.09 0.11 0.09 0.1  0.1  0.09 0.08 0.09 0.13
 0.1  0.1  0.12 0.08 0.11 0.09 0.06 0.13 0.06 0.12 0.07 0.09 0.07 0.13
 0.1  0.07 0.1  0.12 0.21 0.1  0.11 0.09 0.11 0.08 0.08 0.11 0.08 0.11
 0.11 0.19 0.08 0.13 0.11 0.11 0.1  0.16 0.1  0.09 0.06 0.11 0.09 0.07
 0.13 0.05 0.12 0.07 0.06 0.13 0.1  0.15 0.08 0.06 0.1  0.08 0.12 0.08
 0.1  0.1  0.11 0.12 0.11 0.12 0.08 0.08 0.11 0.08 0.14 0.11 0.1  0.1
 0.1  0.07 0.05 0.13 0.14 0.09 0.06 0.13 0.1  0.07 0.18 0.12 0.14 0.12
 0.09 0.13 0.1  0.13 0.08 0.1  0.09 0.1  0.09 0.09 0.08 0.07 0.06 0.11
 0.08 0.11 0.1  0.1  0.07 0.1  0.21 0.13 0.16 0.14 0.1  0.13 0.09 0.15
 0.13 0.07 0.13 0.06 0.05 0.08 0.12 0.11 0.08 0.11 0.08 0.11 0.06 0.12
 0.07 0.11 0.14 0.13 0.15 0.06 0.13 0.11 0.11 0.11 0.15 0.13 0.14 0.05
 0.11 0[...]<jupyter_text>Now we plot the distribution. The easiest first look is a histogram.<jupyter_code>plt.hist(props, bins = 5)
plt.xlabel("proportion")
plt.ylabel("frequency")
plt.show()<jupyter_output><empty_output><jupyter_text>The resolution is a bit inappropriate, given that we deal with integers. To increase the bin number would be a good idea. Maybe we should also plot a confidence interval.<jupyter_code>CI= sm.stats.proportion_confint(n*p, n, alpha=0.05)
print(CI)

plt.axvspan(CI[0],CI[1], alpha=0.2, color='yellow')
plt.hist(props, bins = 50)
plt.xlabel("proportion")
plt.ylabel("frequency")
plt.axvline(p, color="black")<jupyter_output>(0.04120108046379837, 0.15879891953620165)<jupyter_text>In a binomial distribution, the distribution is given by the proportion and the sample size. Therefore we could calculate a confidence interval from one measurement.      How can we now describe the distribution? Summary statistics:<jupyter_code>print("the minimum is:",  min(props))
print("the maximum is:",  max(props))
print(statistics.mean(props))<jupyter_output>0.10157000000000001<jupyter_text>Is the mean a good way to look at our distribution?<jupyter_code>n = 50 # number of trials
p = 0.02 # probability of each trial

s = np.random.binomial(n, p, 1000) #simulation repeating the experiment 1000 times
props = s/n

CI= sm.stats.proportion_confint(n*p, n, alpha=0.05)
print(CI)

plt.axvspan(CI[0],CI[1], alpha=0.2, color='yellow')
plt.hist(props, bins = 20)
plt.xlabel("proportion")
plt.ylabel("frequency")
plt.axvline(p, color="black")
plt.axvline(statistics.mean(props), color="red")

print(statistics.mean(props))<jupyter_output>(0.0, 0.05880530708179099)
0.019880000000000002<jupyter_text>Count data/ the Poisson distribution The Poisson distribution is built on count data, e.g. the numbers of raisins in a Dresdner Christstollen, the number of geese at any given day between Blaues Wunder and Waldschlösschenbrücke, or radioactive decay. So lets use a Geiger counter and count the numbers of decay per min.<jupyter_code>freq =1.6
s = np.random.poisson(freq, 1000)

plt.hist(s, bins = 20)
plt.xlabel("counts per minute")
plt.ylabel("frequency")
plt.axvline(freq, color="black")<jupyter_output><empty_output><jupyter_text>Confidence intervals for a Poisson distribution Similar to the binomial distribution, the distribution is defined by sample size and the mean.  Also for Poisson, one can calculate an also asymmetrical confidence interval:<jupyter_code>freq =1.6
s = np.random.poisson(freq, 1000)
CI = st.poisson.interval(0.95,freq)

plt.axvspan(CI[0],CI[1], alpha=0.2, color='yellow')
plt.hist(s, bins = 20)
plt.xlabel("counts per minute")
plt.ylabel("frequency")
plt.axvline(freq, color="black")<jupyter_output><empty_output><jupyter_text>For a poisson distribution, poisson error can be reduced, when increasing the counting population, in our case lets count for 10 min instead of 1 min, and see what happens.<jupyter_code>CI = np.true_divide(st.poisson.interval(0.95,freq*10),10)
print(CI)
s = np.true_divide(np.random.poisson(freq*10, 1000),10)

plt.axvspan(CI[0],CI[1], alpha=0.2, color='yellow')
plt.hist(s, bins = 70)
plt.xlabel("counts per minute")
plt.ylabel("frequency")
plt.axvline(freq, color="black")<jupyter_output>[0.9 2.4]<jupyter_text>What is the difference between Poisson and Binomial? Aren't they both kind of looking at count data?   Yes, BUT:   Binomial counts events versus another event, e.g. for the cells there are two options, normal versus deformed. A binomial distribution is about comparing the two options.   Poisson counts with an open end, e.g. number of mutations.    Continuous data Let's import the count data you have generated with Robert. When you download it from Google sheets (https://docs.google.com/spreadsheets/d/1Ek-23Soro5XZ3y1kJHpvaTaa1f4n2C7G3WX0qddD-78/editgid=0), it comes with spaces. Try to avoid spaces and special characters in file names, they tend to make trouble.   I renamed it to 'BBBC001.csv'.<jupyter_code>dat = pd.read_csv(os.path.join('data','BBBC001.csv'), header=1, sep=';')

print(dat)<jupyter_output>Annotator name (pseudonym is ok)  BBBC001 manual count  \
0                            Robert                   370   
1                          <NAME>.                   365   
2                          <NAME>.                   390   
3                         <NAME>..                   370   
4                          <NAME>.                   383   
..                              ...                   ...   
83                         <NAME>.                   356   
84                          <NAME>.                   368   
85                         <NAME>.                   367   
86                           Julia                    367   
87                            Marie                   367   

    BBBC001 CLIJ Voronoi Otsu Labeling  BBBC001 StarDist  BBBC001 Find Maxima  
0                                367.0             379.0                  NaN  
1                                360.0             373.0                375.0  
2                          [...]<jupyter_text>For now we will focus on the manual counts, visualise it and perform summary statistics.<jupyter_code>man_count = dat.iloc[:,1].values


plt.hist(man_count,bins=100)<jupyter_output><empty_output><jupyter_text>There are different alternatives of displaying such data, some of which independent of distribution. You will find documentation in the graph galery: https://www.python-graph-gallery.com/<jupyter_code>sns.kdeplot(man_count)<jupyter_output><empty_output><jupyter_text>A density plot is sometimes helpful to see the distribution, but be aware of the smoothing and that you loose the information on sample size.<jupyter_code>sns.stripplot(y=man_count)
sns.swarmplot(y=man_count)
sns.violinplot(y=man_count)<jupyter_output><empty_output><jupyter_text>this plot is useful, but the density function can sometimes be misleading and lead to artefacts dependent on the sample size. Unless explicitely stated, sample sizes are usually normalised and therefore hidden!<jupyter_code>sns.boxplot(y=man_count)<jupyter_output><empty_output><jupyter_text>Be aware that boxplots hide the underlying distribution and the sample size.  So the "safest" plot, when in doubt, is to combine boxplot and jitter:<jupyter_code>ax = sns.swarmplot(y=man_count, color="black")
ax = sns.boxplot(y=man_count,color="lightgrey")<jupyter_output><empty_output><jupyter_text>The boxplot is very useful, because it directly provides non-parametric summary statistics:  Min, Max, Median, Quartiles and therefore the inter-quartile range (IQR). The whiskers are usually the highest point that is within 1.5x the quartile plus the IQR. Everything beyond that is considered an outlier. Whiskers are however not always used in this way! The mean and standard diviation are not visible in a boxplot, because it is only meaningful in distributions that center around the mean. It is however a part of summary statistics:<jupyter_code>dat["BBBC001 manual count"].describe()<jupyter_output><empty_output><jupyter_text>Normal distribution We assume that our distribution is "normal".  First we fit a normal distribution to our data.<jupyter_code>(mu, sigma) = norm.fit(man_count)
n, bins, patches = plt.hist(man_count, 100,density=1)


# add a 'best fit' line
y = norm.pdf(bins, mu, sigma)
l = plt.plot(bins, y, 'r--', linewidth=2)

#plot
plt.xlabel('manual counts counts')
plt.ylabel('binned counts')
plt.title(r'$\mathrm{Histogram\ of\ manual\ counts:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu, sigma))

plt.show()<jupyter_output><empty_output><jupyter_text>Is it really normally distributed? What we see here is already one of the most problematic properties of a normal distribution: The susceptibility to outliers.  In normal distributions the confidence interval is determined by the standard diviation. A confidence level of 95% equals 1.96 x sigma.<jupyter_code>#plot

(mu, sigma) = norm.fit(man_count)
n, bins, patches = plt.hist(man_count, 100,density=1)

# add a 'best fit' line
y = norm.pdf(bins, mu, sigma)
l = plt.plot(bins, y, 'r--', linewidth=2)


plt.xlabel('manual counts counts')
plt.ylabel('binned counts')
plt.title(r'$\mathrm{Histogram\ of\ manual\ counts:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu, sigma))
plt.axvspan((mu-1.96*sigma),(mu+1.96*sigma), alpha=0.2, color='yellow')
plt.axvline(mu, color="black")


plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1000697350069735, 'fraction_numerical': 0.12325662482566248, 'mean_word_length': 2.9947771587743732, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 77, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '13434', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 5087, 'n_tokens_neox': 4018, 'n_words': 1480, 'path': 'biostatistics/stats2.ipynb', 'repo_name': 'Nablaaa/Bio-image_Analysis_with_Python', 'size': '11478'}