// Ce script nécessite le module Atoms CASCI

clear;

// paramètres de la loi de Weibull

beta_forme = 0.845;
eta_echelle = 126;

// génération des données

// 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

// 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 exponentielle

lnR = log(R);

a_exp=sum(tt.*lnR)/sum(tt.^2);

Rexp = 1-cdfexponential(tt, -a_exp);

// tracé

scf(0);
clf;

subplot(2,2,1)
plot(tt, lnR, "o")
xpoly([tt(1), tt($)], [a_exp*tt(1), a_exp*tt($)]);
xtitle("Diagramme semi-logarithmique (loi exponentielle)", "t (j)", "ln R")
xstring(240, -0.2, "$\lambda ="+string(-a_exp)+"$");

// droite de Henry : quantiles loi normale

t_norm = idfnormal(F, 0, 1);
[a_norm, b_norm, sigmanorm] = reglin(tt', t_norm'); // régression linéaire
sigma_norm = 1/a_norm;
mu_norm = -b_norm*sigma_norm;

Rnorm = cdfnormal(tt, mu_norm, sigma_norm);

subplot(2, 2, 2)
plot(tt, t_norm, "o");
xpoly([tt(1), tt($)], [a_norm*tt(1) + b_norm, a_norm*tt($) + b_norm]);
xtitle("Droite de Henry (loi normale)", "t (j)", "quantile")
xstring(10, 1.15, "$\mu ="+string(mu_norm)+"\text{ ; } \sigma ="...
+string(sigma_norm)+"$")

// droite de Henry : quantiles loi log-normale

lnt = log(tt);
[a_lognorm, b_lognorm, sigmalognorm] = reglin(lnt', t_norm');
// régression linéaire
sigma_lognorm = 1/a_lognorm;
mu_lognorm = -b_lognorm*sigma_lognorm;

Rlognorm = 1-cdfnormal(lnt, mu_lognorm, sigma_lognorm);

subplot(2, 2, 3)
plot(lnt, t_norm, "o");
xpoly([lnt(1), lnt($)], [a_lognorm*lnt(1) + b_lognorm, a_lognorm*lnt($) + b_lognorm]);
xtitle("Droite de Henry (loi log-normale)", "ln t", "quantile")
xstring(0.2, 1.15, "$\mu ="+string(mu_lognorm)+"\text{ ; } \sigma ="...
+string(sigma_lognorm)+"$")

// loi de Weibull

Yweib = log(-log(R));
[a_weib, b_weib, sigma_weib] = reglin(lnt', Yweib');
beta_weib = a_weib;
lambda = exp(-b_weib/beta_weib);

Rweib = 1-cdfweibull(tt, beta_weib, lambda);

subplot(2,2,4)
plot(lnt, Yweib, "o")
xpoly([lnt(1), lnt($)], [a_weib*lnt(1) + b_weib, a_weib*lnt($) + b_weib]);
xtitle("Diagramme de Weibull", "t (j)", "ln R")
xstring(0.2, 0.55, "$\beta ="+string(beta_weib)+"\text{ ; } \lambda ="...
+string(lambda)+"$");

scf(1);
clf;

subplot(2,2,1)
plot(tt, R, "o")
plot(tt, Rexp)
xtitle("Loi exponentielle", "t (j)", "R")

subplot(2,2,2)
plot(tt, R, "o")
plot(tt, 1-Rnorm)
xtitle("Loi normale", "t (j)", "R")

subplot(2,2,3)
plot(tt, R, "o")
plot(tt, Rlognorm)
xtitle("Loi log-normale", "t (j)", "R")

subplot(2,2,4)
plot(tt, R, "o")
plot(tt, Rweib)
xtitle("Loi de Weibull", "t (j)", "R")