# import pandas as pan #pour les données :gestion des data.frame
from scipy.stats import bernoulli, binom, geom, poisson
from datetime import datetime
from sympy import isprime
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st
import pandas as pan

###################Exercice 1#########################

# 1
N = 1000
U = st.uniform.rvs(1, 6, size=N)
D = np.floor(U)

print(D)


table = pan.crosstab(D, columns="freq", normalize=True)
print(table)
# 2
plt.bar(table.index, table['freq'], width=0.1)
# Autre solution
table.plot.bar(width=0.1)


# 3
N = 100000
U = st.uniform.rvs(1, 6, size=N)
D = np.floor(U)
table = pan.crosstab(D, columns="freq", normalize=True)
plt.bar(table.index, table['freq'], width=0.1)
# On obtient une meilleure uniformité du dé.
print(table)
# 4


def Dice(f, N):
    U = st.uniform.rvs(1, f, size=N)
    D = np.floor(U)
    return(D)


pan.crosstab(Dice(8, 10000), columns="freq", normalize=True)

###################Exercice 2#########################
# 1


def Restriction(x, k):
    i = 0
    j = 0
    y = x.copy()
    while(j < len(x)):
        while((j < len(x)) and (x[j] > k)):
            j = j+1
        if(j < len(x)):
            y[i] = x[j]
        i = i+1
        j = j+1
    return(y[0:i])

# Autre solution


def Restriction(x, k):
    return(x[x < k+1])


# 2
D = Dice(8, 10000)
R = Restriction(D, 6)
pan.crosstab(R, columns="freq", normalize=True)


# 3
isprime(D)


def prime(x):
    i = 0
    j = 0
    y = x.copy()
    while(j < len(x)):
        while((j < len(x)) and (not (isprime(int(x[j]))))):
            j = j+1
        if(j < len(x)):
            y[i] = x[j]
        i = i+1
        j = j+1
    return(y[0:i])


P = prime(Dice(12, 10000))
pan.crosstab(P, columns="freq", normalize=True)
# On obtient une loi uniforme sur les nombres premiers en dessous de 12
# Cela correspond au fait que conditionnellement à un évènement
# une loi uniforme donne encore une loi uniforme sur l'évènement


###################Exercice 3#########################

n, p = 10, 0.25
ProbasB = binom.pmf(range(n+1), n, p)  # probabilites
# de la distribution binomiale
plt.bar(range(n+1), ProbasB, width=0.05)
plt.show()

# 1
n, p = 10, 0.5
plt.bar(range(n+1), binom.pmf(range(n+1), n, p), width=0.05)
plt.title("Diagramme en bâton de la loi B(10,0.5)")


n, p = 100, 0.25
plt.bar(range(51), binom.pmf(range(51), n, p), width=0.2)
plt.title("Diagramme en bâton de la loi B(100,0.25)")
# Observer la forme de courbe en cloche que l'on expliquera plus tard en cours et au TP4

n = 10
plt.bar(range(n+1), poisson.pmf(range(n+1), 2), width=0.05)
plt.title("Diagramme en bâton de la loi P(2)")


n = 20
plt.bar(range(n+1), poisson.pmf(range(n+1), 10), width=0.2)
plt.title("Diagramme en bâton de la loi P(10)")
# Observer la forme de courbe en cloche que l'on expliquera plus tard en cours et au TP4


n = 5
plt.ylim((0, 1))
plt.bar(range(n), geom.pmf(range(n), 0.75), width=0.05)
plt.title("Diagramme en bâton de la loi G(0.75)")


n = 20
plt.bar(range(n), geom.pmf(range(n), 0.25), width=0.1)
plt.title("Diagramme en bâton de la loi G(0.25)")

# 2 partie simulation

n, p = 10, 0.5
b = binom.rvs(n, p, size=1000)
print(b)
tab = pan.crosstab(b, columns="freq", normalize=True)
print(tab)
plt.bar(tab.index, tab["freq"], width=0.05)
x = tab.index
plt.plot(x, binom.pmf(x, n, p), 'r+')
plt.title("Diagramme en bâton d'une simulation de la loi B(10,0.5)")


mu = 10
p = poisson.rvs(mu, size=100000)
tab = pan.crosstab(p, columns="freq", normalize=True).values
x = range(int(min(tab)), int(min(tab))+len(tab))
plt.bar(x, np.reshape(tab, len(tab)), width=0.1)
plt.plot(x, poisson.pmf(x, mu), 'r+')
plt.title("Diagramme en bâton d'une simulation de la loi P(10)")
plt.show()


# 3

n = 25
x = np.arange(n+1)
plt.bar(x+.6, binom.pmf(x, 50, 0.2), width=0.3, color="green")
plt.bar(x+.3, binom.pmf(x, 250, 0.04), width=0.3, color="blue")
plt.bar(x, binom.pmf(x, 500, 0.02), width=0.3, color="red")
plt.plot(x, poisson.pmf(x, 10), 'o', color="black")
plt.title("Convergence vers la loi de Poisson P(10)")

###################Exercice 4#########################

# 1
n, p = 10, 0.25
plt.step(range(n+1), binom.cdf(range(n+1), n, p))
plt.title("Fonction de répartition de la loi B(10,0.25)")


# 2
n, p = 100, 0.5
plt.step(range(n+1), binom.cdf(range(n+1), n, p))
plt.title("Fonction cumulative de la loi B(100,0.5)")


n = 15
plt.step(range(n+1), poisson.cdf(range(n+1), 3))
plt.title("Fonction cumulative de la loi P(3)")

binom.cdf(3, 50, 0.2)  # P(X\leq 3)=0.005656
1-binom.cdf(30, 50, 0.2)  # P(X> 30)=1.1024325896613618e-10

###################Exercice 5#########################
# 1


def r(x, y):
    z = np.sqrt(x**2+y**2)
    return(z)


r(3, 4)
# on a défini la norme euclidienne (longueur),
# puis on a calculé celle du vecteur (3,4)


# 2
def g(N, p, k):
    nf = np.math.factorial(N)
    kf = np.math.factorial(k)
    nkf = np.math.factorial(N-k)
    z = nf/(kf*nkf)*(p**k)*((1-p)**(N-k))
    return(z)


g(20, 0.5, 10)
st.binom.pmf(10, n=20, p=0.5)  # [1] 0.1761971
x = range(21)
plt.bar(x, height=[g(20, 0.5, k) for k in x], width=0.1)
plt.plot(x, st.binom.pmf(x, n=20, p=0.5), '.', color="red")

###################Exercice 6#########################
# 1
N = 10**4
t0 = datetime.now().microsecond
s = 0
for i in range(1, N):
    if (i % 2 == 0):
        s = s+i**2
# Pour tous les entiers pairs (congrus à 0 modulo 2),
# on ajoute leur carré, donc on fait la somme des
# carrés des entiers pairs entre 1 et N
t1 = datetime.now().microsecond
t = np.arange(1, N)
s2 = np.sum(t[t % 2 == 0]**2)
t2 = datetime.now().microsecond
print("Valeurs obtenues :\n s = ", s, " temps : ",
      t1-t0, "\n s2 = ", s2, "temps : ", t2-t1, "\n")
# Réponse chez moi:
# Valeurs obtenues :
#  s =  166716670000  temps :   10595 microsecondes
# s2 =  166716670000 temps :  1186 microsecondes

# On a comparé la vitesse de calcul de la somme
# des carrés d'entiers pairs entre 1 et 10000,
# soit par une boucle, soit en utilisant les
# opérations sur les vecteurs de numpy, et la deuxième
# est 10 fois plus rapide...


# 2 et 3
N = 2*10**5  # 2k+1,k dans 0:99 ou  2k-1, k dans 1:100
t0 = datetime.now().microsecond
s = 0
for i in range(1, N):
    if (i % 2 == 1):
        s = s+1/i**2
t1 = datetime.now().microsecond
t = np.arange(100000)
s2 = np.sum(1/(2*t+1)**2)
t2 = datetime.now().microsecond
print("Valeurs obtenues :\n s = ", s, " temps : ",
      t1-t0, "\n s2 = ", s2, "temps : ", t2-t1, "\n")
# Valeurs obtenues :
# s =  1.2336980501361898  temps :  0:00:00.045075
# s2 =  1.2336980501361696 temps :  0:00:00.000537

# La deuxième méthode est 100 fois plus rapide !
# 3 Bilan: il vaut mieux éviter les boucles en python!!