Décision dans l'incertain¶

pi, wi, W, title = choisir_loi(poisson) # poisson, binomiale, discrete

fig = plt.gcf()

#plt.subplot(1,3,1);
plt.bar(wi,pi,width=1)
plt.title("loi"); 
plt.show()

#plt.subplot(1,3,2);
plt.step(wi,np.cumsum(pi),where='post',marker='o',markerfacecolor='r');
plt.title("cdf :" + title)
plt.show()

#plt.subplot(1,3,3);
hist_plot(W,wi)
plt.title("loi empirique"); 
plt.show()

c=10;cm=20;cs=5;cf=200;

# la fonction coût j(u,w)
def jj(u,w):
    # A vous : écrire la fonction j (voir slides)


def J(u,pi,wi):
    # A vous : écrire la fonction J = la moyenne de j (voir slides)


def resultats(pi,wi,dessin='oui'):
# draw the function and the minimum 

    # remplir val avec les valeurs de J(u) quand u varie de -10 a 100 avec un pas de 1
    u=np.linspace(-10,100,num=111)
    val=np.zeros(u.size);
    for i in range(0,u.size):
        # A vous : val[i]= ??
        # val[i]=  # modifier
    
    # recherche du min 
    # A vous : trouver la valeur uopt
    # uopt = ?
    
    if dessin=='oui':
        plt.plot(u,val);
        plt.title("fonction coût");
        plt.plot(uopt,val[imin],marker='o',markerfacecolor='r')

        print("Nombre optimal de journaux a commander *%f*, " % uopt,end='');
        print("Moyenne de la demande *%f*\n" % sum(wi*pi));

    return uopt


pi, wi, W, title = choisir_loi(discrete) # poisson, binomiale, discrete
print('Cas discrêt: ',end='')
resultats(pi,wi)
plt.show()
    
pi, wi, W, title = choisir_loi(binomiale) # poisson, binomiale, discrete
print('Cas binomial: ',end='')
resultats(pi,wi)
plt.show()

pi, wi, W, title = choisir_loi(poisson) # poisson, binomiale, discrete
print('Cas Poisson: ',end='')
resultats(pi,wi)
plt.show()

pi, wi, W, title = choisir_loi(binomiale) # poisson, binomiale, discrete

# On dessine F 
fstar = (cm-c)/(cm+cs) # valeur limite
Fv= np.cumsum(pi) # calcul de la fonction de repartition
plt.clf()
xx=np.hstack((0,wi)) # on rajoute 0 a wi
yy=np.hstack((0,Fv)) # on rajoute 0 a Fv
plt.step(xx,yy,where='post') # tracé de la fonction de répartition
plt.plot(xx,fstar*np.ones(xx.size)) # droite horizontale de valeur fstar

# chercher l'indice qui donne uopt avec la formule du cours
kopt=0
for i in range(0,Fv.size):
    if Fv[i] >= fstar:
        kopt=i
        break;

plt.step(xx,yy,where='post') # marker='o',markerfacecolor='r');
plt.plot(xx,fstar*np.ones(xx.size))
plt.plot(wi[kopt],Fv[kopt],marker='o',markeredgewidth=5,markeredgecolor='r',markerfacecolor='r');

wi=np.array([10,50,80])
pi=np.array([3/4,1/8,1/8])

resultats(pi,wi)

plt.clf()

def draw_cost(u,W):
    #count, bins, ignored =   plt.hist(j(u,W),density=True)
    samples=jj(u,W)
    max_samples=max(samples)
    hist_plot(samples)
    return max_samples

# poisson , discrete, binomiale
pi, wi, W, title = choisir_loi(binomiale)

u=49
# commence par traiter le problème par simulation (avec W)
max_samples=draw_cost(u,W);
plt.show()

# calcul de la loi exacte
valeurs=jj(u,wi) # toutes les valeurs de j(u,wi[i])
support = np.sort(list(dict.fromkeys(valeurs))) 
      # on calcule le support de la loi
      # en enlevant les valeurs en double
loi=np.zeros(support.size) # on va calculer la loi de j(u,W)
for i in range(support.size):
    for k in range(wi.size):
        if jj(u,wi[k]) == support[i]:
            loi[i]=loi[i]+pi[k]

# on tronque l'histogramme au delà de la valeur max_samples
for imax in range(support.size):
    if support[imax] > max_samples:
        print(imax);break
imax=imax-1
            
# On trace cet histogramme
plt.bar(support[0:imax], loi[0:imax],width=5)

# On regarde maintenant un cas ou le vendeur à déjà des journaux 
# et ou il paye un cout fixe s'il commande des journaux 
# on voit une stratégie [s,S] 

def Jtilde(u,x,pi,wi):
    return cf*(u>0) + J(u+x,pi,wi) -c*x

# Le but est de verifier que l'on a une stratégie [s,S]

# Coix du modèle : poisson , discrete, binomiale
pi, wi, W, title = choisir_loi(poisson)

# Calcul de S
# S=uopt, où uopt est l'argmin de J(u) (cf Slide 11/20)
# uopt est caclulé dans la fonction résultat
uopt=resultats(pi,wi,dessin='non') # on ne veux pas de dessin ici
S=uopt 

# Il nous faut maintenant calculer s
xv=np.linspace(0,2*max(wi), num=200+1);
xuopt=xv.copy();
U=np.linspace(0,2*max(wi), num=200+1); # 
Ju=U.copy();

# A vous : pour chaque valeur du stock xv[i], calculer le cout avec la commande U[j]
# puis chercher la commande optimale que vous rangez dans xuopt[i]
# for i in range(0, xv.size):
#    for j in range(0, U.size):

# s est la valeur la plus grande où le contrôle est non nul
# on cherche donc l'indice iopt de la valeur la plus petite où le 
# contrôle est non nul. xv[iopt-1] donne la valeur de s cherchée.
# A vous : 
#iopt=-1;
#for i in range(0, xuopt.size):
# (iopt-1) est le dernier indice ou la controle est non nul
# et s est la valeur de x correspondante xv[iopt-1]
s = xv[iopt-1]
print("1) Calcul par une méthode brutale : valeur de s=%f et de S=%f\n" % (s, S))

plt.clf()

# on vérifie que pour x en dessous de s, x+uopt=cte=S
plt.plot(xv[:int(s)],xv[:int(s)]+xuopt[:int(s)])

# au dela xuopt=0 (on ne commande rien), donc on est sur la droite y = x!
plt.plot(xv[int(s)+1:150],xv[int(s)+1:150]+xuopt[int(s)+1:150])

plt.show()

# On vérifie le formule de s donnée dans le cours (p.11 cours-4) fonctionne
xv= np.linspace(0,2*max(wi), num=200+1);
Jv=np.zeros(xv.size);
for i in range(0,xv.size):
    Jv[i]=J(xv[i],pi,wi)
JS=J(S,pi,wi) # J(S)
costs = Jv - (cf + JS);
# calcul du plus grand z où J(z) >= c_F + J(S)
iopt=-1
for i in range(0,xv.size):
    if costs[i] <=0:
        iopt=i;
        break;
if iopt < 0: 
    s=0 
else:
    s=xv[iopt-1]
    
print("2) Formule du cours: valeur de s=%f et de S=%f\n" % (s, S))