Kumaraswamy distribution

Author

Affiliation

Beniamino Sartini

University of Bologna

library(dplyr)
# required for figures  
library(ggplot2)
library(gridExtra)
# required for render latex 
library(backports)
library(latex2exp)
# Random seed 
set.seed(1)

In probability and statistics, the Kumaraswamy’s double bounded distribution is a family of continuous probability distributions defined on the interval \((0,1)\). It is similar to the Beta distribution, but much simpler to use especially in simulation studies since its probability density function, cumulative distribution function and quantile functions can be expressed in closed form. This distribution was originally proposed by Poondi Kumaraswamy. The distribution function is in Equation 1, the density function is in Equation 2, the quantile in Equation 3 and finally the random generator.

1 Distribution function

For \(x \in (0,1)\) and \(a > 0\), \(b > 0\) the distribution function is defined as:

\[ F(x; a, b) = 1 - (1-x^a)^{b} \tag{1}\]

In R code:

\[ \begin{aligned} \text{pkumaraswamy}({} & x, \; a = 1, \; b = 1, \\ & \text{log.p} = {\color{red}{\text{FALSE}}}, \\ & \text{lower.tail} = {\color{green}{\text{TRUE}}} \end{aligned} \]

grid <- seq(0, 1, 0.01)
bind_rows(
  tibble(label = "a = 0.5, b = 0.5", a = 0.5, b = 0.5, x = grid, p = pkumaraswamy(x, a, b)),
  tibble(label = "a = 5, b = 1", a = 5, b = 1, x = grid, p = pkumaraswamy(x, a, b)),
  tibble(label = "a = 1, b = 3", a = 1, b = 3, x = grid, p = pkumaraswamy(x, a, b)),
  tibble(label = "a = 2, b = 2", a = 2, b = 2, x = grid, p = pkumaraswamy(x, a, b)),
  tibble(label = "a = 2, b = 5", a = 2, b = 5, x = grid, p = pkumaraswamy(x, a, b))
) %>%
  ggplot()+
  geom_line(aes(x, p, color = label))+
  theme_bw()+
  labs(x = TeX("$x$"), y = TeX("$F(x;a,b)$"), color = NULL)+
  theme(legend.position = "top")

Figure 1: Kumaraswamy distribution function for different parameters.

R

pkumaraswamy <- function(x, a = 1, b = 1, log.p = FALSE, lower.tail = TRUE){
  x[x<=0] <- 0
  x[x>=1] <- 1
  p <- 1 - (1 - x^a)^b
  if (!lower.tail){
    p <- 1 - p
  }
  if (log.p){
    p <- base::log(p)
  }
  return(p)
}

Python

import numpy as np
def pkumaraswamy(x, a=1, b=1, log_p=False, lower_tail=True):
    x = np.where(x >= 1, 0, np.where(x <= 0, 1, x))
    p = 1 - (1 - x**a)**b
    if not lower_tail:
      p = 1 - p
    if log_p:
      p = np.log(p)
    return p

2 Density function

For \(x \in (0,1)\) and \(a > 0\), \(b > 0\) the density function is defined as:

\[ f(x; a, b) = abx^{a-1} (1-x^a)^{b-1} \tag{2}\]

In R code:

\[ \begin{aligned} \text{dkumaraswamy}({} & x, \; a = 1, \; b = 1, \\ & \text{log.p} = {\color{red}{\text{FALSE}}}) \end{aligned} \]

grid <- seq(0.01, 0.99, 0.01)
bind_rows(
  tibble(label = "a = 0.5, b = 0.5", a = 0.5, b = 0.5, x = grid, p = dkumaraswamy(x, a, b)),
  tibble(label = "a = 5, b = 1", a = 5, b = 1, x = grid, p = dkumaraswamy(x, a, b)),
  tibble(label = "a = 1, b = 3", a = 1, b = 3, x = grid, p = dkumaraswamy(x, a, b)),
  tibble(label = "a = 2, b = 2", a = 2, b = 2, x = grid, p = dkumaraswamy(x, a, b)),
  tibble(label = "a = 2, b = 5", a = 2, b = 5, x = grid, p = dkumaraswamy(x, a, b))
) %>%
  ggplot()+
  geom_line(aes(x, p, color = label))+
  theme_bw()+
  labs(x = TeX("$x$"), y = TeX("$f(x;a,b)$"), color = NULL)+
  theme(legend.position = "top")

Figure 2: Kumaraswamy density function for different parameters.

R

dkumaraswamy <- function(x, a = 1, b = 1, log.p = FALSE){
  p <- a*b*x^(a - 1)*(1 - x^a)^(b - 1)
  p[x<0|x>1] <- 0
  if (log.p){
    p <- base::log(p)
  }
  return(p)
}

Python

import numpy as np
def dkumaraswamy(x, a=1, b=1, log_p=False):
    p = a * b * (x ** (a - 1)) * ((1 - x ** a) ** (b - 1))
    p = np.where((x < 0) | (x > 1), 0, p)
    if log_p:
        p = np.log(p)
    return p

3 Quantile function

For \(p \in (0,1)\) and \(a > 0\), \(b > 0\) the quantile function is defined as:

\[ F^{-1}(p; a, b) = 1 - (1-p^{\frac{1}{b}})^{\frac{1}{a}} \tag{3}\]

In R code:

\[ \begin{aligned} \text{qkumaraswamy}({} & p, \; a = 1, \; b = 1, \\ & \text{log.p} = {\color{red}{\text{FALSE}}}, \\ & \text{lower.tail} = {\color{green}{\text{TRUE}}} \end{aligned} \]

R

qkumaraswamy <- function(p, a = 1, b = 1, log.p = FALSE, lower.tail = TRUE){
  if (log.p) {
    p <- exp(p)
  }
  if (!lower.tail){
    p <- 1 - p
  }
  x <- (1 - (1-p)^(1/b))^(1/a)
  return(x)
}

Python

import numpy as np
def qkumaraswamy(p, a = 1, b = 1, log_p = False, lower_tail = False):
    p = np.where(p >= 1, 1, np.where(p <= 0, 0, p))
    if log_p: 
       p = np.exp(p)
    if not lower_tail: 
       p = 1 - p
    x = (1 - (1 - p)**(1/b))**(1/a)
    return x

4 Random generator

In general, given a random variable \(X \sim F_{X}\) following a certain distribution \(F_{X}\) and the correspondent quantile function \(F_{X}^{-1}\) (i.e. \(F_{X}^{-1}(F(X)) = X\)) it is possible to generate a random sample by simulating the grades. In steps:

• Step 1.: simulate the grades \(u \sim \text{Unif}[0,1]\).

• Step 2.: apply the quantile function of the desired distribution on the grades to obtain a simulations.

In R code:

\[ \text{rkumaraswamy}(n, \; a = 1, \; b = 1) \]

params <- c(a = 1.5, b = 1.5)
x <- rkumaraswamy(5000, a = params[1], b = params[2])

ker <- density(x, from = min(x), to = max(x))
ker$pdf_th <- dkumaraswamy(ker$x, a = params[1], b = params[2])
#ker$pdf_th <- ker$pdf_th/sum(ker$pdf_th)
ker$pdf_emp <- ker$y#/sum(ker$y)

ggplot()+
  geom_line(aes(ker$x, ker$pdf_th), color = "red")+
  geom_line(aes(ker$x, ker$pdf_emp))+
  theme_bw()+
  labs(x = TeX("$x$"), y = TeX("$f(x;a,b)$"), color = NULL)+
  theme(legend.position = "top")

Figure 3: Kumaraswamy density function for different parameters.

R

rkumaraswamy <- function(n, a = 1, b = 1){
  u <- runif(n, min = 0, max = 1)
  x <- qkumaraswamy(u, a = a, b = b, log.p = FALSE, lower.tail = TRUE)
  return(x)
}

Python

import numpy as np
def rkumaraswamy(n, a = 1, b = 1): 
    u = np.random.uniform(0, 1, n)
    x = qkumaraswamy(u, a = a, b = b, log_p = False, lower_tail = True)
    return x