nate — Jun 24, 2014, 4:47 PM
# ASM 13th edition, 2.3
par(mfrow=c(1,1))
#### Terminology ####
# d: corresponds to a pdf (density function)
# p: corresponds to a cdf (distribution function OR probability)
# c: corresponds to a normalizing constant
#### Parameters ####
uniform.params = list(d=0, u=1)
beta.params = list(theta=1,a=1,b=0.5)
exponential.params = list(theta=5)
weibull.params = list(t=1, theta=1)
gamma.parms = list(a=0.5, theta=100)
pareto.params = list(a=0.5, theta=5)
lognormal.params = list(mu=2, sigma=1.9555)
#### Uniform Distribution ####
# d<x<u, u>d>=0
duniform <- function(x, d=uniform.params$d, u=uniform.params$u) {
x^0* (1/(u-d))
}
puniform <- function(x, d=uniform.params$d, u=uniform.params$u) {
if(x<=d) {
0
} else if(d<=x && x<=u) {
(x-d)/(u-d)
} else { # x>=u
1
}
}
Euniform <- function(d=uniform.params$d, u=uniform.params$u) {
(d+u)/2
}
Vuniform <- function(d=uniform.params$d, u=uniform.params$u) {
(u-d)^2/12
}
# recognize uniform distribution by finite support and lack of x in the density function
#### Beta Distribution ####
# 0<x<theta, a>0, b>0, theta>0
# theta=a=b=1 is a uniform distribution
dbeta <- function(x, theta=beta.params$theta, a=beta.params$a, b=beta.params$b) {
c = gamma(a+b)/(gamma(a)*gamma(b))
c*x^(a-1)*(theta-x)^(b-1)
}
Ebeta <- function(theta=beta.params$theta, a=beta.params$a, b=beta.params$b) {
theta*a/(a+b)
}
Vbeta <- function(theta=beta.params$theta, a=beta.params$a, b=beta.params$b) {
theta^2*(a*b)/((a+b)^2*(a+b+1))
}
# mode = theta*(a-1)/(a+b-2) if a>1 && b>1
# distribution function can be evaluated if a or b is an integer
# note that distribution function requires numerical integration (incomplete beta function)
#### Exponential Distribution ####
# x>=0, theta>0
dexponential <- function(x, theta=exponential.params$theta) {
c=1/theta
c*exp(-x/theta)
}
pexponential <- function(x, theta=exponential.params$theta) {
1 - exp(-x/theta)
}
Eexponential <- function(theta=exponential.params$theta) {
theta
}
Vexponential <- function(theta=exponential.params$theta) {
theta^5
}
# the higher the parameter, the more weight placed on higher numbers
#### Weibull Distribution ####
# x>=0, t>0, theta>0
dweibull <- function(x, t=weibull.params$t, theta=weibull.params$theta) {
c=t*theta^(-t)
c*x^(t-1)*exp(-(x/theta)^t)
}
# evaluating the moments requires evaluating the gamma function (numerical techniques)
Eweibull <- function(t=weibull.params$t, theta=weibull.params$theta) {
theta*gamma(1+1/t)
}
# E(X^2)
E2Weibull <- function(t=weibull.params$t, theta=weibull.params$theta) {
theta*gamma(1+2/t)
}
# the distribution has a non-zero mode when t>1.
# notice that the distribution with t=0.5 puts a lot of weight on small numbers
# to make up for this, it will also put higher weight that very large numbers
#### Gamma Distribution ####
# x>=0, a>0, theta>0
dgamma <- function(x, a=gamma.params$a, theta=gamma.params$theta) {
c = 1/(gamma(a)*theta^a)
c*x^(a-1)*exp(-x/theta)
}
Egamma <- function(x, a=gamma.params$a, theta=gamma.params$theta) {
a*theta
}
Vgamma <- function(x, a=gamma.params$a, theta=gamma.params$theta) {
a*theta^2
}
# when a is an integer, a gamma random variable with parameters a and theta is
# a sum of a independent random variables with parameter theta. when a = 1, the
# gamma random variable is exponential. the gamma distribution is called an erlang
# distribution when a is an integer
#
# you recognize a gamma distribution when the density function has e raised to a
# multiple of x and in addition has x raised to a power. contrast this to weibull
# where e is raised to a multiple of a power of x
#
# the distribution function may be evaluated if a is an integer; otherwise numerical
# techniques are needed
#### Pareto Distribution ####
# x>=0, theta>0, a>0
dpareto <- function(x,a=pareto.params$a, theta=pareto.params$theta) {
c=a*theta^a
c/(theta+x)^(a+1)
}
Epareto <- function(a=pareto.params$a, theta=pareto.params$theta) {
theta/(a-1) # a>1 or moment doesnt exist
}
# E(X^2)
E2pareto <- function(a=pareto.params$a, theta=pareto.params$theta) {
2*theta^2/((a-1)*(a-2)) # a>2 or moment doesnt exist
}
# you recognize a Pareto when the density function has a denominator with x plus
# a constant raised to a power
#
# the mean can be infinite for pareto distribution
#### Single-Parameter Pareto Distribution ####
# x>=theta, a>0
# theta is not considered a parameter since it must be selected in advance
# based on what you want the range to be
dsinglepareto <- function(x,a=pareto.params$a, theta=pareto.params$theta) {
c=a*theta^a
c/x^(a+1)
}
Esinglepareto <- function(a=pareto.params$a, theta=pareto.params$theta) {
a*theta/(a-1) # a>1 or moment doesnt exist
}
# E(X^2)
E2singlepareto <- function(a=pareto.params$a, theta=pareto.params$theta) {
a*theta^2/(a-2) # a>2 or moment doesnt exist
}
# you recognize a single-parameter Pareto by having support not starting at 0,
# and by the density function have a denominator with x raised to a power
# a beta distribution may also have x raised to a negative power, but it would have
# finite support
#
# a single-parameter Pareto X is a two-parameter Pareto Y shifted by theta: X=Y+theta
# thus it has the same variance, and the mean is theta greater than the mean of a
# two-parameter pareto with the same parameters
#### Lognormal distribution ####
dlognormal <- function(x, mu=lognormal.params$mu, sigma=lognormal.params$sigma) {
z=(log(x)-mu)/sigma
exp(-z^2/2)/(x*sigma*sqrt(2*pi))
}
Elognormal <- function(mu=lognormal.params$mu, sigma=lognormal.params$sigma) {
exp(mu+0.5*sigma^2)
}
E2lognormal <- function(mu=lognormal.params$mu, sigma=lognormal.params$sigma) {
exp(2*mu + 2*sigma^2)
}
# you recognize a lognormal distribution by the ln x in the exponent
# if Y is normal, then X = e^y is lognormal with the same parameters u and sigma
# thus, to calculate the distribution function use
# dnorm(z) where z = (log(x)-u)/sigma and dnorm is the cdf of the standard normal distribution
# more generally E(X^k) = E(e^(k*Y)) = My(k) where My(k) is the moment generating
# function of the corresponding normal distribution
#
# the mode is exp(u-sigma^2)
# as sigma gets lower, the distribution flattens out
#### Plots ####
# uniform
uniform.params = list(d=0, u=10)
plot(duniform, uniform.params$d, uniform.params$u, col='red', lwd=1, main="Uniform distribution", ylab="f(x)")
legend('topleft', c('d=0, u=1'),
lty=1, lwd=1, col=c('red', 'blue'), bty='n', cex=.9)
# beta
beta.params = list(theta=1, a=1, b=0.5)
plot(dbeta, 0, 1, col='red', lwd=1, main="Beta distribution", ylab="f(x)")
beta.params = list(theta=1, a=2, b=1)
plot(dbeta, 0, 1, add=TRUE, col='blue', lwd=1)
beta.params = list(theta=1, a=6, b=3)
plot(dbeta, 0, 1, add=TRUE, col='green', lwd=1)
beta.params = list(theta=1, a=18, b=9)
plot(dbeta,0,1, add=TRUE, col='brown', lwd=1)
legend('topleft', c('a=1, b=0.5','a=2, b=1','a=6, b=3','a=18, b=9') ,
lty=1, lwd=1, col=c('red', 'blue', 'green',' brown'), bty='n', cex=.9)
# exponential
exponential.params = list(theta=25)
plot(dexponential, 0, 100, col='red', lwd=1, main="Exponential distribution", ylab="f(x)")
exponential.params = list(theta=50)
plot(dexponential, 0, 100, add=TRUE, col='blue', lwd=1)
exponential.params = list(theta=100)
plot(dexponential, 0, 100, add=TRUE, col='green', lwd=1)
legend('topright', c('theta=25', 'theta=50', 'theta=100'),
lty=1, lwd=1, col=c('red','blue', 'green'), bty='n', cex=.9)
# weibull
weibull.params = list(t=0.5, theta=25)
plot(dweibull, 0, 100, col='red', lwd=1, main="Weibull distribution", ylab="f(x)")
weibull.params = list(t=1,theta=50)
plot(dweibull, 0, 100, add=TRUE, col='blue', lwd=1)
weibull.params = list(t=2, theta=100/sqrt(pi))
plot(dweibull, 0, 100, add=TRUE, col='green', lwd=1)
legend('topright', c('t=0,5 theta=25', 't=1, theta=50', 't=2, theta=100/sqrt(pi)'),
lty=1, lwd=1, col=c('red','blue', 'green'), bty='n', cex=.9)
# gamma
gamma.params = list(a=0.5, theta=100)
plot(dgamma, 0, 100, col='red', lwd=1, main="Gamma distribution", ylab="f(x)")
gamma.params = list(a=5,theta=10)
plot(dgamma, 0, 100, add=TRUE, col='blue', lwd=1)
gamma.params = list(a=50, theta=1)
plot(dgamma, 0, 100, add=TRUE, col='green', lwd=1)
legend('topright', c('a=0.5, theta=100','a=5, theta=10','a=50, theta=1') ,
lty=1, lwd=1, col=c('red', 'blue', 'green'), bty='n', cex=.9)
# pareto
pareto.params = list(a=0.5, theta=5)
plot(dpareto, 0, 100, col='red', lwd=1, main="Pareto distribution", ylab="f(x)")
pareto.params = list(a=2, theta=50)
plot(dpareto, 0, 100, add=TRUE, col='blue', lwd=1)
pareto.params = list(a=5, theta=200)
plot(dpareto, 0, 100, add=TRUE, col='green', lwd=1)
legend('topright', c('a=0.5, theta=5','a=2, theta=50','a=5, theta=200') ,
lty=1, lwd=1, col=c('red', 'blue', 'green'), bty='n', cex=.9)
# lognormal
lognormal.params = list(mu=2, sigma=1.9555)
plot(dlognormal, 0, 100, col='red', lwd=1, main="Lognormal distribution", ylab="f(x)")
lognormal.params = list(mu=3, sigma=1.3506)
plot(dlognormal, 0, 100, add=TRUE, col='blue', lwd=1)
lognormal.params = list(mu=3.5, sigma=0.9078)
plot(dlognormal, 0, 100, add=TRUE, col='green', lwd=1)
legend('topright', c('mu=2, sigma=1.9555', 'mu=3, sigma=1.3506', 'mu=3.5, sigma=0.9078'),
lty=1, lwd=1, col=c('red','blue','green'), bty='n', cex=.9)
