//*********************************************
// BRÛLAGE
// Texte Scopos vol 11, p. 175.
//*********************************************

// On considère que T suit une loi de Weibull avec gamma=0 et  eta=1.
// Choix de beta : le taux de défaillance vaut r(t)= beta t^{beta-1},
// le brûlage n'a de sens que si r(t) est décroissante, il faut donc beta<1.
beta=0.6; // avec gamma=0 et  eta=1 (pas utilisés).

// Autres paramètres :
alpha=0.8;
n=20;      // taille de l'échantillon pour la méthode de la fdr empirique
N=100;     // nombre d'échantillons
k=10000;   // nombre de points pour calcul de l'intégrale Phi(u)

////////////////////////////////////////////////////////////////////////
// Durée optimale du brûlage : cas F connue
// Tracé de la courbe y=Phi(u) pour u dans [0,1] :
//   F^{-1}(u) = (-log(1-u)) .^ (1/beta)
//   Phi(u) = 1/mu * \int_0^F^{-1}(u) R(s) ds
mu=gamma(1 + 1/beta); // espérance de la loi Weibull(beta)
// Le domaine d'intégration pour le calcul de Phi est [0, +oo[,
// on calcule les t=F^{-1}(u) pour u dans [0,1] (k points)
u=[1:k-1] / k;
t=(-log(1-u)) .^ (1/beta);
// Il est facile de calculer une valeur approchée de l'intégrale 
// Psi(t)= \int_0^u R(s) ds
// par la méthode de Riemann (rectangles) :
R= exp(- t .^beta);   // valeurs de la fonction
dt= t - [0 t(1:k-2)]; // largeur des rectangles
rect= R .* dt;        // aires des rectangles
phi= cumsum(rect);
// N.B. : comme R(t) décroît, phi minore l'intégrale et
Phi= cumsum([1 R(1:k-2)] .* dt);
// majore l'intégrale.
// On a donc une majoration de l'erreur sur les valeurs de Phi
// (augmenter k si l'erreur est trop importante).
// La valeur approchée de l'intégrale par la méthode des trapèzes est
// (Psi+psi)/2, celle de Phi est (Psi+psi)/2mu.
printf('\n');
printf('Majorant erreur sur Phi : %f', max(Phi-phi)/mu);
Phi=(Phi+phi)/2/mu;
// Tracé Phi en rouge (style=5) Phi(0)=0, Phi(1)=1.
xbasc();
xselect();
plot2d([0 u 1], [0 Phi 1], style=5*ones(t));
// Remarque : Phi est convexe pour beta < 1, concave pour beta > 1
// et Phi(u) = u si beta = 1 (évident sur la déf. de Phi).
//
// Calcul du b optimal :
[m,I]=max((alpha - Phi) ./ (ones(u)-u));
// u^*=u(I) et  b^*=F(u^*)=1-R(I)
printf('\n');
printf('Durée optimale théorique brûlage  : b = %.2f (u =%.2f)', t(I), u(I));
// On trace la sécante optimale et l'abcisse u^*
// sur le graphe en marron (style=25)
xsegs([1      1-R(I)], [alpha Phi(I)], 25);
xsegs([1-R(I) 1-R(I)], [ 0    Phi(I)], 25);
////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////
// Durée optimale du brûlage : cas F estimée par la f.d.r. empirique
// de l'échantillon.

// On simule N n-échantillons de loi de Weibull avec beta<1 
// (pour taux de panne r décroissant) : R(t)=1-F(t)=exp(-t^beta ).

// Tirage des échantillons : F continue strict. croissante, d'où
// si U suit une loi uniforme, F^{-1}(U) suit une loi de Weibull(beta).
U=zeros(N,n);
U=rand(U,'uniform');      // N n-échantillons uniformes
T= (-log(U)) .^ (1/beta); // N n-échantillons Weibull(beta)
// On réordonne chaque échantillon par ordre croissant :
T= -sort(-T,'c');

// On note v au lieu de u, et Phi_n(v) la fonction empirique
// (Phi_n(v) est constante par morceaux).
v=[1:n] / n;
// Matrice des accroissements : T_1, T_2-T_1, ... , T_n -T_{n-1} :
Phi_n = T - [zeros(N,1) T(:,[1:n-1])];
// Matrice n*T_1, (n-1)*(T_2-T_1), ... , 1*(T_n -T_{n-1}) :
Phi_n = Phi_n .* (ones(N,1) * [n:-1:1]);
// On somme :
Phi_n = cumsum(Phi_n, 'c');
// On divise chaque ligne par la somme des T_i :
Phi_n = Phi_n ./ (sum(T,'c') * ones(1,n));
// Tracé de la fonction Phi_n (en bleu : 2) pour le 1er échantillon 
// (1ère ligne de la matrice).
w=[0:n-1] / n; // Décaler les abcisses pour plot2d2
plot2d2(w, Phi_n(1,:), style=2*ones(t), strf="000");
// Matrice des pentes pour les N échantillons : 
// alpha/1, (alpha-Phi_n(1/n))/(1-1/n), ... (alpha-Phi_n((n-1)/n)/(1-(n-1)/n))
Num=[alpha*ones(N,1), [alpha*ones(N,n-1) - Phi_n(:,[1:n-1])]];
Den= ones(N,1) * (ones(w) - w);
Pentes = Num ./ Den;
[M,J]=max(Pentes, 'c'); 
// J est compris entre 1 et n, l'indice optimal entre 0 et n-1,
// on "corrige" J : J = J-1
J=J-ones(J);
// u et b "optimaux" : u^*==v(J)   b^*=T_{J}
// On trace la sécante optimale et l'abcisse u^* pour le 1er échantillon
// en vert (style=13).
xsegs([1       v(J(1))], [alpha Phi_n(1,J(1))],13);
xsegs([v(J(1)) v(J(1))], [ 0    Phi_n(1,J(1))],13);
//
printf('\n');
printf('Durée optimale brûlage (1er éch.) : b = %.2f\n',T(1,J(1)));
// Moyenne et écart-type empiriques pour les durées estimées de brûlage
// sur les N échantillons :
B=zeros(N,1);
for j=1:N
  if J(j)==0 then
    B(j)=0;
  else
    B(j)=T(j,J(j));
  end
end;
Bmoy= sum(B) / N;
VarB= sum((B-Bmoy).^2) / (N-1);
printf('Moyenne sur %d échantillons      : b = %.2f\n', N, Bmoy);
printf('Écart-type sur %d échantillons   : e = %.2f\n', N, sqrt(VarB));