• variabile casuale di Kumaraswamy e Beta

Kumaraswamy

pkuma <- function(x,a,b){     return(1-(1-x^a)^b)}
rkuma <- function(n,a,b){return((1-(1-runif(n))^(1/b))^(1/a))}
qkuma <- function(q,a,b){return((1-(1-q)^(1/b))^(1/a))}

dkuma <- function(x,a,b){     return(a*b*x^(a-1)*(1-x^a)^(b-1))}
modakuma <- function(a,b){return(max(0,min(1,((a-1)/(a*b-1))^(1/a))))}

dx <- 0.001;
x <- seq(0,1,dx)
a <- 4
b <- 1.5
m<-b*gamma(1+1/a)*gamma(b)/gamma(1+1/a+b)
m2<-b*gamma(1+2/a)*gamma(b)/gamma(1+2/a+b)
s2<-m2-m^2

p <- m*(m*(1-m)/s2-1)
q <- p*(1-m)/m



png("c:/temp/Grafico_Cumulate_Kumaraswamy_vs_Beta.png",width=800,heigh=600)
plot(x,pkuma(x,a,b),
     type='l',
     axes=FALSE,xlab=NA,ylab=NA,
     ylim=c(0,1), xlim=c(0,1),
     main="Kumaraswamy vs. Beta\nFunzioni di ripartizione",
     col=1
     );
abline(0,0);lines(c(0,0),c(-10,100))
abline(1,0,lty=2)
lines(x,pbeta(x,p,q),col=2)
legend(0.145,1,col=i,
       legend=c(paste("v.c. di Kumaraswamy\na =",a,"  \nb =",b),paste("\nv.c. Beta\np =",p,"  \nq =",q)),
       lty=1,bty='n',
       ncol=1)
text(0,1,"1",adj=1.4)
text(1,0,"1",adj=c(0.5,1.4))
dev.off()