R Markdown

Severity distributions

We can start with the exponential distribution [

f_X(x)= e^{-x} ]

x<-c(0:50)
x<-c(0:50)
plot(x,dexp(x,0.005),xlab="x",ylab="f(x)",type="l",col="blue")

plot(x,dexp(x,0.01),xlab="x",ylab="f(x)",type="l",col="blue")

plot(x,dexp(x,0.05),xlab="x",ylab="f(x)",type="l",col="blue")

plot(x,dexp(x,0.1),xlab="x",ylab="f(x)",type="l",col="blue")

plot(x,dexp(x,1),xlab="x",ylab="f(x)",type="l",col="blue")

The Gamma distribution

Remember that this distribution has the following density function:

\[ f_X(x)= \frac{( x/ \theta)^\alpha}{x \Gamma( \alpha)} e^{-x / \theta}, \quad x>0, \alpha >0, \theta>0 \]

scaleparam <- seq(100, 250, by = 50)
shapeparam <- 2:5
x <- seq(0, 1000, by = 1)
fgamma <- dgamma(x, shape = 2, scale = scaleparam[1])
plot(x, fgamma, type = "l", ylab = "Gamma Density")
for(k in 2:length(scaleparam)){
  fgamma <- dgamma(x,shape = 2, scale = scaleparam[k])
  lines(x,fgamma, col = k)
}
legend("topright", c("scale=100", "scale=150", "scale=200", "scale=250"), lty=1, col = 1:4)

# Varying Shape Gamma Densities
fgamma <- dgamma(x, shape = shapeparam[1], scale = 100)
plot(x, fgamma, type = "l", ylab = "Gamma Density")
for(k in 2:length(shapeparam)){
  fgamma <- dgamma(x,shape = shapeparam[k], scale = 100)
  lines(x,fgamma, col = k)
}
legend("topright", c("shape=2", "shape=3", "shape=4", "shape=5"), lty=1, col = 1:4)

Remember that we can move the shape to \(\alpha=1\)? We can check if that happens:

scaleparam <- seq(100, 250, by = 50)
x <- seq(0, 1000, by = 1)
fgamma <- dgamma(x, shape = 1, scale = scaleparam[1])
plot(x, fgamma, type = "l", ylab = "Gamma Density")
for(k in 2:length(scaleparam)){
  fgamma <- dgamma(x,shape = 1, scale = scaleparam[k])
  lines(x,fgamma, col = k)
}
legend("topright", c("scale=100", "scale=150", "scale=200", "scale=250"), lty=1, col = 1:4)

For the example using Gamma(3,0.5)

x<-c(0:25)
plot(x,dgamma(x,4,0.5),xlab="x",ylab="f(x)",type="l",col="blue")

pgamma(4,4,0.5)-pgamma(2,4,0.5)
## [1] 0.1238884

For the second example

# Using normal distribution:
1-pnorm(2.191)
## [1] 0.0142259

Weibull distribution.

par(mfrow=c(1, 2), mar = c(4, 4, .1, .1))

# Varying Scale Weibull Densities
z<- seq(0,400,by=1)
scaleparam <- seq(50,200,50)
shapeparam <- seq(1.5,3,0.5)
plot(z, dweibull(z, shape = 3, scale = scaleparam[1]), type = "l", ylab = "Weibull density")
for(k in 2:length(scaleparam)){
  lines(z,dweibull(z,shape = 3, scale = scaleparam[k]), col = k)}
legend("topright", c("scale=50", "scale=100", "scale=150", "scale=200"), lty=1, col = 1:4)

# Varying Shape Weibull Densities
plot(z, dweibull(z, shape = shapeparam[1], scale = 100), ylim=c(0,0.012), type = "l", ylab = "Weibull density")
for(k in 2:length(shapeparam)){
  lines(z,dweibull(z,shape = shapeparam[k], scale = 100), col = k)}
legend("topright", c("shape=1.5", "shape=2", "shape=2.5", "shape=3"), lty=1, col = 1:4)

1-pweibull(600,shape=0.5,scale=300)
## [1] 0.2431167
pweibull(500,0.5,300)
## [1] 0.7250028

Pareto distribution

require(actuar)
## Loading required package: actuar
## Warning: package 'actuar' was built under R version 4.1.2
## 
## Attaching package: 'actuar'
## The following objects are masked from 'package:stats':
## 
##     sd, var
## The following object is masked from 'package:grDevices':
## 
##     cm
par(mfrow=c(1, 2), mar = c(4, 4, .1, .1))

# Varying Scale Pareto Densities
x <- seq(1, 3000, by = 1)
scaleparam <- seq(2000, 3500, 500)
shapeparam <- 1:4

# varying the shape parameter
plot(x, actuar::dpareto(x, shape = shapeparam[1], scale = 2000), ylim=c(0,0.002),type = "l", ylab = "Pareto density")
for(k in 2:length(shapeparam)){
  lines(x, actuar::dpareto(x, shape = shapeparam[k], scale = 2000), col = k)
}
legend("topright", c(expression(alpha~'=1'), expression(alpha~'=2'), expression(alpha~'=3'), expression(alpha~'=4')), lty=1, col = 1:4)

# Varying Shape Pareto Densities
plot(x, actuar::dpareto(x, shape = 3, scale = scaleparam[1]), type = "l", ylab = "Pareto density")
for(k in 2:length(scaleparam)){
  lines(x, actuar::dpareto(x, shape = 3, scale = scaleparam[k]), col = k)
}
legend("topright", c(expression(theta~'=2000'), expression(theta~'=2500'), expression(theta~'=3000'), expression(theta~'=3500')), lty=1, col = 1:4)

For the example with Pareto

a<-1-actuar::ppareto(75, shape = 3, scale = 100)
a
## [1] 0.1865889
b<-1-actuar::ppareto(50, shape = 3, scale = 100)
a/b
## [1] 0.6297376