Simple Birnbaum Saunders
# BS
A <- function(t, alpha, beta){
(t+beta)/(2*alpha*sqrt(t^3)*sqrt(beta))
}
a <- function(t, alpha, beta){
(1/alpha)*(sqrt(t/beta)-sqrt(beta/t))
}
phi1 <- function(t) {
dnorm(a(t,alpha, beta),0,1)*A(t,alpha,beta)
}
FBS <- function(t){
pnorm(a(t,alpha, beta),0,1)
}
# jpeg(filename = 'myplot1.jpg', units = 'cm', width = 24, height = 14, res=800)
par(mfrow=c(1,2))
# change alpha
# pdf
alpha=0.5
beta=1
curve(phi1, from=0, to=6, col=2, lwd=3, lty=2, ylim=c(0,1),xlab='t',ylab='PDF')
alpha=1
beta=1
curve(phi1, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
alpha=1.5
beta=1
curve(phi1, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
legend('topright', c(expression(alpha==0.5~','~beta==1),expression(alpha==1.0~','~beta==1),expression(alpha==1.5~','~beta==1)),
col=2:4, lty=2:4, lwd=3, bty='n')
# cdf
alpha=0.5
beta=1
curve(FBS, from=0, to=6, col=2, lwd=3, lty=2, ylim=c(0,1),xlab='t',ylab='CDF')
alpha=1
beta=1
curve(FBS, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
alpha=1.5
beta=1
curve(FBS, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
legend('bottomright', c(expression(alpha==0.5~','~beta==1),expression(alpha==1.0~','~beta==1),expression(alpha==1.5~','~beta==1)),
col=2:4, lty=2:4, lwd=3, bty='n')
par(mfrow=c(1,2))
# change beta
# pdf
alpha=0.5
beta=1
curve(phi1, from=0, to=6, col=2, lwd=3, lty=2, ylim=c(0,1),xlab='t',ylab='PDF')
alpha=1
beta=2
curve(phi1, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
alpha=1.5
beta=3
curve(phi1, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
legend('topright', c(expression(alpha==0.5~','~beta==1),expression(alpha==0.5~','~beta==2),expression(alpha==0.5~','~beta==3)),
col=2:4, lty=2:4, lwd=3, bty='n')
# cdf
alpha=0.5
beta=1
curve(FBS, from=0, to=6, col=2, lwd=3, lty=2, ylim=c(0,1),xlab='t',ylab='CDF')
alpha=1
beta=2
curve(FBS, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
alpha=1.5
beta=3
curve(FBS, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
legend('bottomright', c(expression(alpha==0.5~','~beta==1),expression(alpha==0.5~','~beta==2),expression(alpha==0.5~','~beta==3)),
col=2:4, lty=2:4, lwd=3, bty='n')
Finite Mixture Birnbaum Saunders
# FM - Birnbaum-Saunders
A <- function(t, alpha, beta){
(t+beta)/(2*alpha*sqrt(t^3)*sqrt(beta))
}
a <- function(t, alpha, beta){
(1/alpha)*(sqrt(t/beta)-sqrt(beta/t))
}
phi1 <- function(t,alpha,beta) {
dnorm(a(t,alpha, beta),0,1)*A(t,alpha,beta)
}
f1 <- function(t){
if(sum(pii)!=1) stop('The constraints have not met!')
v <- 0
for (i in 1:k) {
alpha=alphaa[i]
beta=betaa[i]
v <- v + pii[i]*phi1(t,alpha=alphaa[i],beta=betaa[i])
}
v
}
skew-t-normal Finite Mixture Birnbaum Saunders
g <- function(x, nu){
s1=gamma((nu+1)/2)/(sqrt(nu*pi)*gamma(nu/2))
s2=(1 + (x^2)/nu)^(-(nu+1)/2)
}
psi <- function(t){
2*g(a(t,alpha,beta), nu)*pnorm(lambda*a(t,alpha,beta))*A(t,alpha,beta)
}
### TBS
g <- function(t, nu){
dt(x=t, df=nu)
}
phi2 <- function(t,alpha,beta,nu){
g(a(t=t,alpha=alpha,beta=beta), nu=nu)*A(t=t,alpha=alpha,beta=beta)
}
f2 <- function(t){
if(sum(pii)!=1) stop('The constraints have not met!')
v <- 0
for (i in 1:k) {
alpha = alphaa[i]
beta = betaa[i]
nu = nuu[i]
v <- v + pii[i]*phi2(t=t,alpha=alphaa[i],beta=betaa[i],nu=nuu[i])
}
v
}
### SNBS
phi3 <- function(t,alpha,beta,lambda){
2*dnorm(a(t,alpha,beta))*pnorm(lambda*a(t,alpha,beta))*A(t,alpha,beta)
}
f3 <- function(t){
if(sum(pii)!=1) stop('The constraints have not met!')
v <- 0
for (i in 1:k) {
alpha <- alphaa[i]
beta <- betaa[i]
lambda <- lambdaa[i]
v <- v + pii[i]*phi3(t,alpha,beta,lambda)
}
v
}
### STNBS
g <- function(x, nu){
s1=gamma((nu+1)/2)/(sqrt(nu*pi)*gamma(nu/2))
s2=(1 + (x^2)/nu)^(-(nu+1)/2)
}
psi <- function(t,alpha,beta,nu,lambda){
2*g(a(t,alpha,beta), nu)*pnorm(lambda*a(t,alpha,beta))*A(t,alpha,beta)
}
f4 <- function(t){
if(sum(pii)!=1) stop('The constraints have not met!')
v <- 0
for (i in 1:k) {
alpha <- alphaa[i]
beta <- betaa[i]
nu <- nuu[i]
lambda <- lambdaa[i]
v <- v + pii[i]*psi(t,alpha=alphaa[i],beta=betaa[i],nu=nuu[i],lambda=lambdaa[i])
}
v
}
Visualization of The Distribution
par(mfrow=c(2,2))
alphaa=c(0.5,0.75)
betaa=c(1,2)
pii=c(0.3,0.7)
k=2
lambdaa=c(1,3)
nuu=c(2,3)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(lambda==(1~','~3)),expression(nu==(2~','~3))))
#----------------------------------------------------------
lambdaa=c(-1,1)
nuu=c(1,1)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(lambda==(-1~','~1)),expression(nu==(1~','~1))))
#----------------------------------------------------------
lambdaa=c(-2,2)
nuu=c(2,3)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(lambda==(-2~','~2)),expression(nu==(2~','~3))))
#----------------------------------------------------------
lambdaa=c(-4,4)
nuu=c(4,6)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(lambda==(-4~','~4)),expression(nu==(4~','~6))))
# dev.off()
Visualization of The Distribution
par(mfrow=c(2,2))
alphaa=c(0.5,0.75)
betaa=c(1,2)
pii=c(0.3,0.7)
lambdaa=c(-4,4)
nuu=c(4,6)
k=2
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(alpha==(0.5~','~0.75)),expression(beta==(1~','~2))))
#----------------------------------------------------------
alphaa=c(1,0.75)
betaa=c(1,2)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(alpha==(1~','~0.75)),expression(beta==(1~','~2))))
#----------------------------------------------------------
alphaa=c(1,2)
betaa=c(1,3)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(alpha==(1~','~2)),expression(beta==(1~','~3))))
#----------------------------------------------------------
alphaa=c(2,1)
betaa=c(3,1)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('topleft', c(expression(alpha==(2~','~1)),expression(beta==(3~','~1))))
Visualization of The Distribution
par(mfrow=c(2,2))
alphaa=c(0.5,0.75)
betaa=c(1,2)
pii=c(0.3,0.7)
lambdaa=c(-4,4)
nuu=c(4,6)
k=2
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('right', c(expression(pi==(0.3~','~0.7))))
#----------------------------------------------------------
pii=c(0.4,0.6)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('right', c(expression(pi==(0.4~','~0.6))))
#----------------------------------------------------------
pii=c(0.6,0.4)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.9),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('right', c(expression(pi==(0.6~','~0.4))))
#----------------------------------------------------------
pii=c(0.7,0.3)
curve(f1, from = 0.1, to=6, col=2, lwd=3, lty=2, ylim=c(0,1.7),xlab='t',ylab='PDF')
curve(f2, from = 0.1, to=6, col=3, lwd=3, lty=3, add=T)
curve(f3, from = 0.1, to=6, col=4, lwd=3, lty=4, add=T)
curve(f4, from = 0.1, to=6, col=5, lwd=3, lty=5, add=T)
legend('topright', c('FM-BS','FM-TBS','FM-SNBS','FM-STNBS'),
col=2:5, lty=2:5, lwd=3, bty='n')
legend('right', c(expression(pi==(0.5~','~0.75))))
Visualization of The Distribution
fY = function(t){
p*phi1(t, alpha1, beta1)+(1-p)*phi1(t, alpha2, beta2)
}
alpha1=0.5; alpha2=1
beta1=1; beta2=2
p = 0.6
curve(fY, from=0, to=6, col=2, lwd=3, lty=2, ylim=c(0,1),xlab='t',ylab='Density')
alpha1=0.6; alpha2=5
beta1=2; beta2=0.5
p = 0.6
curve(fY, from=0, to=6, col=2, lwd=3, lty=2, ylim=c(0,1),xlab='t',ylab='Density')
follow me on RPubs for more!