Chapter3, GFGMO figures
n=10; p=1; q=1; k=5; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.07, c("Blomq", "rho", "tau"), cex=0.5, col=c("red", "blue",
"green"), lty=1:3)
##############################
n=10; p=1; q=2; k=1; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.65, -0.07, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
#####################################
n=10; p=1; q=3; k=5; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.004, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
##################################
n=10; p=2; q=1; k=1; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.3, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
#############################
n=10; p=2; q=2; k=2; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.11, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
##############################
n=10; p=2; q=3; k=1; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.05, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
#############################
n=10; p=3; q=1; k=1; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.45, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
###############################
n=10; p=3; q=2; k=5; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.15, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
#############################
n=10; p=3; q=3; k=5; b=k/n; d=(n-k)/n; betta=0.1
#Blomq
measure=function(theta) (theta*(b*betta + d))*( 1- ((1/2)^p) ) ^ (2*q)
plot(measure, type="l", col="red",xlab=expression(theta),xlim=c(-1,1), lty=1)
#rho
measure=function(theta) (((12/(p^2))*theta*(b*betta + d))
) * ((beta(2/p, q+1))^2)
plot(measure, type="l", col="blue", xlim=c(-1,1), main="betta=0.1",lty=1,
add=TRUE)
#tau
measure=function(theta) (4/(p^2)) * ( (theta*(b*betta + d)
)*((beta(2/p, q+1))^2)
+ (theta*(b*betta + d))*(( beta(2/p, q) - (1+(p*q))*beta((2/p) +1, q) )^2)
+ ((theta*(b*betta + d))^2)*( beta(2/p, 2*q) - (1+(p*q)
)*beta((2/p) +1, 2*q) )^2 )
plot(measure, type="l", col="green",xlim=c(-1,1), main="betta=0.1",lty=3,
add=TRUE)
legend(0.5, -0.15, c("Blomq", "rho", "tau"), cex=0.5,
col=c("red", "blue", "green"), lty=1:3)
Chapter3, MGFGMO figures
ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*((1-(v^p))^q)*
(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1; n=10; p=1; q=1.5
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*(log(1+theta*(1-2*u)*
(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+ (log(1+(theta*beta)*(
(1-((u)^p))^(q-1))*((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
############################
#plots
betta=0.1; k=1; n=10; p=1; q=1.5
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p,q) - (1+(p*q))
*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2
+ 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),
lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n))
plot(measure, type="l", col="blue",xlim=c(-1,1), lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)* ( (1/(p^2))*
(beta((4/p)-1, 2*q +1))+((1/p)*(beta(3/p, 2*q +1)))) )+
(4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1), lty=3,
add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q)) +
(theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1), lty=4,add=TRUE)
legend(0.5, -0.1, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=2; n=10; p=1; q=1.5
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2+(4/(p^2))*theta*betta*(k/n)*( beta(2/p,q) -
(1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) - (1+p*q)*
beta(2/p+1, 2*q) )^2
+ 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1), lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)* ( (1/(p^2))*
(beta((4/p)-1, 2*q +1)) + ((1/p)*(beta(3/p, 2*q +1)))) ) +
(4/15*theta*(n-k)/n))
plot(measure, type="l", col="yellow",xlim=c(-1,1), lty=3,
add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q))
+ (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1), lty=4,add=TRUE)
legend(0.5, -0.1, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=1; q=1.5
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p,q) - (1+(p*q))*
beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) - (1+p*q)*
beta(2/p+1, 2*q) )^2
+ 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),
lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1), lty=2,
add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)* ( (1/(p^2))*
(beta((4/p)-1, 2*q +1)) + ((1/p)*(beta(3/p, 2*q +1)))) ) +
(4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1), lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q)) +
(theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1), lty=4,add=TRUE)
legend(0.5, -0.05, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*((1-(v^p))^q)*
(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1; n=10; p=1; q=2
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#plots
betta=0.1; k=1; n=10; p=1; q=2
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q))
+ (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.09, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=1; q=2
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q))
+ (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"),
cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
@
<<>>=
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*((1-(v^p))^q)*
(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1; n=10; p=1; q=3
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m)
{sl1=sl1+ (log(1+(theta*beta)*
((1-((u)^p))^(q-1))*((1-((v)^p))^(q-1))*(1-(1+p*q)*
((u)^p))*(1-(1+p*q)*((v)^p))))*
(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
ظ«plots
betta=0.1; k=2; n=10; p=1; q=3
#tau
measure=function(theta) (4/(p^2))*theta*betta*
(k/n)*(beta(2/p, q+1))^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p,q) - (1+(p*q))*
beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) - (1+p*q)*
beta(2/p+1, 2*q) )^2 +
2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1), lty=2,
add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)* ( (1/(p^2))*
(beta((4/p)-1, 2*q +1)) + ((1/p)*(beta(3/p, 2*q +1)))) ) +
(4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1), lty=3,
add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q)) +
(theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1), lty=4,add=TRUE)
legend(0.5, -0.09, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=1; q=3
#tau
measure=function(theta) (4/(p^2))*theta*betta*
(k/n)*(beta(2/p, q+1))^2 +(4/(p^2))*theta*betta*
(k/n)*( beta(2/p,q) - (1+(p*q))*
beta((2/p)+1,q) )^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p, 2*q) -(1+p*q)*beta(2/p+1, 2*q) )^2
+2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) +
((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1), lty=2,
add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)* ( (1/(p^2))*
(beta((4/p)-1, 2*q +1)) + ((1/p)*(beta(3/p, 2*q +1)))) ) +
(4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1), lty=3,
add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q)) +
(theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1), lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*((1-(v^p))^q)*
(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1;n=10
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#plots
betta=0.1; k=2; n=10; p=1.5; q=1
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=1.5; q=1
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*
((1-(v^p))^q)*(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1;n=10
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#plots
betta=0.1; k=1; n=10; p=1.5; q=2
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.1, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=1.5; q=2
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*((1-(v^p))^q)*
(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1; n=10
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#####
measures of dependence
betta=0.1; k=1; n=10; p=1.5; q=3
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#ginni
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
B=theta*beta*(k/n)*(1-(1/2)^p)^(2*q) + theta/4*(n-k)/n
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.1, c("tau", "rho", "gamma", "Blomq"), cex=0.5,
col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=1.5; q=3
#tau
measure=function(theta) (4/(p^2))*theta*betta*
(k/n)*(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*
(k/n)*( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1),lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#ginni
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"),
lty=1:4)
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*
((1-(v^p))^q)*(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999, beta=beta,
theta=theta)$root
}
return(x)
}
k=1; n=10; p=2; q=1
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#plots
betta=0.1; k=1; n=10; p=2; q=1
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*(k/n)*
((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*((1-(1/2)^p)^(2*q)) +
(theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.09, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=2; q=1
#tau
measure=function(theta) (4/(p^2))*theta*betta*
(k/n)*(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*
(k/n)*( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) +
((1/p)*(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.07, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"),
lty=1:4)
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*
((1-(v^p))^q)*(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999,
beta=beta, theta=theta)$root
}
return(x)
}
k=1; n=10; p=2; q=2
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*
((u)^p))*(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#plots
betta=0.1; k=1; n=10; p=2; q=2
#tau
measure=function(theta) (4/(p^2))*theta*betta*(k/n)*
(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*(k/n)*
( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1), lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.07, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
betta=0.1; k=5; n=10; p=2; q=2
#tau
measure=function(theta) (4/(p^2))*theta*betta*
(k/n)*(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*
(k/n)*( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.07, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"), lty=1:4)
> ff=function(theta,beta,k,n,p,q) {
x=matrix(0,n,2)
u=c();vprin=c()
for(i in 1:n){
u[i]=runif(1);vprin[i]=runif(1)
f=function(v,theta, beta) {theta*((n-k)/n)*v*(1-v)*(1-2*u[i])+
v + (k/n)*v*theta*beta*((1-(u[i])^p)^(q-1))*
((1-(v^p))^q)*(1- ((u[i])^p*(1+p*q)))} - vprin[i]
f=Vectorize(f)
x[i,1]=u[i]
x[i,2]=uniroot(f,lower=0.00001, upper=0.99999,
beta=beta, theta=theta)$root
}
return(x)
}
k=1; n=10; p=2; q=3
l=function(theta, beta) - sum({
A1=0
B1=0
for (ii in 0:n) {A1=A1+ choose(n,ii)*
(log(1+theta*(1-2*u)*(1-2*v)))^(n-ii)
sl1=0
m=ii-n+1
if (m<=0) B1==0
else {for (l in 0:m) {sl1=sl1+
(log(1+(theta*beta)*((1-((u)^p))^(q-1))*
((1-((v)^p))^(q-1))*(1-(1+p*q)*((u)^p))*
(1-(1+p*q)*((v)^p))))*(l/factorial(l))}
B1=B1+ sl1*factorial(ii-n+1)}
}
A1*B1
})
#plots
betta=0.1; k=5; n=10; p=2; q=3
#tau
measure=function(theta) (4/(p^2))*theta*betta*
(k/n)*(beta(2/p, q+1))^2 + (4/(p^2))*theta*betta*
(k/n)*( beta(2/p,q) - (1+(p*q))*beta((2/p)+1,q) )^2 +
(4/(p^2))*theta*betta*(k/n)*( beta(2/p, 2*q) -
(1+p*q)*beta(2/p+1, 2*q) )^2 + 2/9*theta*(n-k)/n
plot(measure, type="l", col="green",xlab=expression(theta),
xlim=c(-1,1), lty=1)
#rho
measure=function(theta) ( (12/(p^2))*theta*betta*
(k/n)*((beta(2/p , q+1))^2) + ((theta/3)*(n-k)/n) )
plot(measure, type="l", col="blue",xlim=c(-1,1),
lty=2, add=TRUE)
#gamma
measure=function(theta) ( (4*theta*betta*(k/n)*
( (1/(p^2))*(beta((4/p)-1, 2*q +1)) + ((1/p)*
(beta(3/p, 2*q +1)))) ) + (4/15*theta*(n-k)/n) )
plot(measure, type="l", col="yellow",xlim=c(-1,1),
lty=3, add=TRUE)
#Blomq
measure=function(theta) ( theta*betta*(k/n)*
((1-(1/2)^p)^(2*q)) + (theta/4*(n-k)/n) )
plot(measure, type="l", col="red",xlim=c(-1,1),
lty=4,add=TRUE)
legend(0.5, -0.06, c("tau", "rho", "gamma", "Blomq"),
cex=0.5, col=c("green", "blue", "yellow", "red"),
lty=1:4)