// Ce script nécessite le module Atoms CASCI

clear;

/// **********
// génération des données
// **********

// Données originales : paramètres de la loi de Weibull

// beta_forme = 0.845;
// eta_echelle = 126;

// Y = linspace(0, 1, 30)
// Y = Y(2:$-1);
// t_orig = floor(idfweibull(Y, beta_forme, eta_echelle))';

t = [2, 5, 9, 13, 17, 22, 27, 39, 39, 39, 52, 64, 64, 76, 86, 97, 108, 121,...
   135, 151, 169, 191, 215, 245, 282, 332]';
t_complet = [t ; 365 ; 365];
N = 28
Ncomplet = size(t_complet, "*");

// nombre cumulé

i = 1;
j = 1
n0 = 1;
nt = size(t, "*");
while i < nt
    if t(i)<>t(i+1) then
        dn(j) = n0; n0 = 1;
        tt(j) = t(i);
        j = j + 1;
    else
        n0 = n0 + 1;
    end
    i = i+1;
end
dn(j) = n0;
tt(j) = t(i);

ndn = j;

n(1) = dn(1);
for i = 2:ndn
    n(i) = n(i-1) + dn(i);
end

// Fréquences cumulées

F = n/(N+1);
R = 1-F;

// loi de Fisher et intervalles de confiance

indices = 1:ndn;

nu1inf = 2*(Ncomplet - indices +1);
nu2inf = 2*indices;
for i = indices
    fisherinf(i) = idffisher(0.95, nu1inf(i), nu2inf(i));
end
Finf = indices'./(indices' + (Ncomplet - indices' + 1).*fisherinf);

nu1sup = 2*(indices + 1);
nu2sup = 2*(Ncomplet - indices);
for i = indices
    fishersup(i) = idffisher(0.95, nu1sup(i), nu2sup(i));
end
Fsup = ((indices' + 1).*fishersup)...
    ./(indices' + (Ncomplet - indices' + 1).*fishersup);

Rsup = 1-Finf;
Rinf = 1-Fsup;

// taux de defaillance lambda
 
for i = 1:ndn-1
    tauxdef(i) = dn(i+1)/(N - n(i))/(tt(i + 1) - tt(i));
end
 
MTTF = mean(t_complet);

// **********
// Tracé
// **********

scf(0);
clf;

subplot(2, 1, 1)
plot(tt, [F, R])
xtitle("Fréquences cumulées", "t (j)", "Probabilité")
legend("F", "R")

subplot(2, 1, 2)
plot(tt(1:$-1), tauxdef)
xtitle("Taux de défaillance", "t (j)", "$\lambda (\mathrm{j}^{-1})$")

scf(1);
clf;

subplot(2,1,1)
plot(tt, Finf, "--", tt, Fsup, "--", tt, F);
xtitle("Courbe de mortalité", "t (j)", "Fréquence cumulée")
legend("min", "max", "F", 2)
a1 = gca();
a1.data_bounds=[0,0;365,1];
a1.auto_scale="off";
a1.tight_limits="on";

subplot(2,1,2)
plot(tt, Rinf, "--", tt, Rsup, "--", tt, R);
xtitle("Courbe de survie (fiabilité)", "t (j)", "Fréquence cumulée")
legend("min", "max", "R")
a2 = gca();
a2.data_bounds=[0,0;365,1];
a2.auto_scale="off";
a2.tight_limits="on";