Bathtub curve and Weibull distribution

bathtub_simulation.m

clear all; close all;

f = @(m,eta,g,t) m/eta^m*(t-g).^(m-1);

t=0.1:0.1:2;
m=0.5; eta=3; g = 0;
lambda = f(m,eta,g,t);
subplot(311);
plot(t,lambda);
hold on;

t=2:0.1:6;
m=1; eta=5; g = 1;
lambda = f(m,eta,g,t);
plot(t,lambda);

t=6:0.1:10;
m=3.5; eta=7; g = 1;
lambda = f(m,eta,g,t);
plot(t,lambda);
title('lambda');

t=0.1:0.1:2;
m=0.5; eta=3; g = 0;
syms x;
f = @(m,eta,g) m/eta^m*(x-g).^(m-1);
int_lambda(x) = int(f(m,eta,g),0,x);
R(x) = exp(-int_lambda(x));
subplot(312);
plot(t,R(t));
hold on;
last = int_lambda(t(end));
subplot(313);
pdf(x)=-diff(R);
plot(t,pdf(t));hold on;

t=2:0.1:6;
m=1; eta=5; g = 1;
int_lambda(x) = last + int(f(m,eta,g),t(1),x);
R(x) = exp(-int_lambda(x));
subplot(312);
plot(t,R(t));
last = int_lambda(t(end));
subplot(313);
pdf(x)=-diff(R);
plot(t,pdf(t));

t=6:0.1:10;
m=3.5; eta=7; g = 1;
int_lambda(x) = last + int(f(m,eta,g),t(1),x);
R(x) = exp(-int_lambda(x));
subplot(312);
plot(t,R(t));
subplot(313);
pdf(x)=-diff(R);
plot(t,pdf(t));

subplot(311);title('Failure rate lambda');
subplot(312);title('Reliability R(t)');
subplot(313);title('Failure PDF f(t)');
sgtitle('Bathtub curve')