6.2 - pdf method for monotone transformations

The pdf method for monotone transformations

The pdf method follows from a specific application of the CDF method
Suppose \(Y\) has pdf \(f_Y(y)\)
Define transformation \(U = g(Y)\), with \(g(Y)\) monotone increasing or decreasing over the support of \(Y\)
Then:

\[f_U(u) = f_Y\left(g^{-1}(u)\right) \left| \frac{d}{du} g^{-1}(u)\right|\]

Monotone functions: definition

A function \(g(Y)\) is monotone increasing over the support of \(Y\) if, for any \(y_2>y_1\) with both \(f_Y(y_1)\), \(f_Y(y_2) > 0\):

\[g(y_2) > g(y_1)\]

A function \(g(Y)\) is monotone decreasing over the support of \(Y\) if, for any \(y_2>y_1\) with both \(f_Y(y_1)\), \(f_Y(y_2) > 0\):

\[g(y_2) < g(y_1)\]

Monotone increasing example

Suppose \(Y\sim BETA\)
Consider \(U = g(Y) = Y^2\):

Proof of the pdf method for monotone increasing transformations

Let \(U= g(Y)\) be a monotone increasing function over the support of \(Y\)
Proceed with the CDF method:

\[\scriptsize F_U(u) = P(U\le u) = P(g(Y) \le u) = P\left(Y \le g^{-1}(u)\right) = F_Y\left(g^{-1}(u)\right)\]

\[ \scriptsize \Rightarrow f_U(u) = \frac{d}{du} F_U(u) = \frac{d}{du} F_Y\left(g^{-1}(u)\right)= f_Y\left(g^{-1}(u)\right)\cdot \frac{d}{du} g^{-1}(u)\]

\[ \scriptsize = f_Y\left(g^{-1}(u)\right)\cdot \left|\frac{d}{du} g^{-1}(u)\right|\]

Monotone decreasing example

Suppose \(Y\sim BETA\)
Consider \(U = g(Y) = e^{-Y}\):

Proof of the pdf method for monotone increasing transformations

Let \(U= g(Y)\) be a monotone decreasing function over the support of \(Y\)
Proceed with the CDF method:

\[\scriptsize F_U(u) = P(U\le u) = P(g(Y) \le u) = P\left(Y \ge g^{-1}(u)\right)\]

\[\scriptsize = 1-F_Y\left(g^{-1}(u)\right)\]

\[ \scriptsize \Rightarrow f_U(u) = \frac{d}{du} F_U(u) = \frac{d}{du}\left(1- F_Y\left(g^{-1}(u)\right)\right)\]

\[\scriptsize = -f_Y\left(g^{-1}(u)\right)\cdot \frac{d}{du} g^{-1}(u)\]

\[ \scriptsize = f_Y\left(g^{-1}(u)\right)\cdot \left|\frac{d}{du} g^{-1}(u)\right|\]

Example: squaring a beta

Suppose \(Y\sim BETA(2,1)\), i.e.:

\[f_Y(y) = \begin{cases} 2y & 0 < y < 1 \\ 0 & otherwise \end{cases}\]

Let \(U = g(Y)= Y^2\). Find \(f_U(u)\).
Support: \(0 < U < 1\)

\[Y=g^{-1}(U) = \sqrt{U} \Rightarrow \frac{d}{du}g^{-1}(u) = \frac{1}{2\sqrt{u}}\]

\[f_U(u) = f_Y(\sqrt{u})\cdot \frac{1}{2\sqrt{u}} = 2\sqrt{u} \cdot \frac{1}{2\sqrt{u}} = \begin{cases} 1 & 0 < u < 1 \\ 0 & otherwise \end{cases} \]

\[\Rightarrow U \sim UNIF(0,1)!\]

Simulating squared betas

library(tidyverse)

(sim_df <- data.frame(Y = rbeta(10000, 2, 1)) 
          %>% mutate(U = Y^2) 
          %>% mutate(fY = dbeta(Y, 2, 1),
                 fU = dunif(U, 0, 1))
) %>% head

          Y         U        fY fU
1 0.7837054 0.6141942 1.5674108  1
2 0.6167265 0.3803516 1.2334531  1
3 0.6942055 0.4819213 1.3884110  1
4 0.8371793 0.7008692 1.6743586  1
5 0.5920146 0.3504812 1.1840291  1
6 0.3868317 0.1496388 0.7736634  1

Plotting results

library(patchwork) # To add plots side-by-side!

ggplot(data = sim_df, aes(x = Y)) + 
    geom_histogram(aes(y=after_stat(density)), 
                   fill = 'goldenrod',color ='black',
                   center = 0.02, binwidth = 0.04)+
    geom_line(aes(y = fY), linewidth = 1) + 
    xlim(c(0,1)) + ylim(c(0,2.5))  + 
    labs(x = 'y', y = expression(f[Y](y))) + 
    theme_classic(base_size = 14)  +


ggplot(data = sim_df, aes(x = U)) + 
    geom_histogram(aes(y=after_stat(density)), 
                   fill = 'cornflowerblue',color ='black',
                   center = -0.98,  binwidth = 0.04)+
    geom_line(aes(y = fU), 
                  linewidth = 1) + 
    xlim(c(0,1)) + ylim(c(0,2.5))  + 
    labs(x = 'u', y = expression(f[U](u))) + 
    theme_classic(base_size = 14)

Simulated densities of Y and U with analytic densities superimposed

Example: exponentiating a beta

Suppose \(Y\sim BETA(2,1)\), i.e.:

\[f_Y(y) = \begin{cases} 2y & 0 < y < 1 \\ 0 & otherwise \end{cases}\]

Let \(U = g(Y)= e^{-Y}\). Find \(f_U(u)\).
Support: \(e^{-1}< U < 1\)

\[Y=g^{-1}(U) = -\ln(U) \Rightarrow \frac{d}{du}g^{-1}(u) = -\frac{1}{u}\]

\[f_U(u) = f_Y(-\ln(u))\cdot \left|-\frac{1}{u}\right| = \begin{cases} \frac{-2\ln(u)}{u} & e^{-1} < u < 1 \\ 0 & otherwise \end{cases} \]

Simulating exponentiated betas

(sim_df <- data.frame(Y = rbeta(10000, 2, 1)) 
          %>% mutate(U = exp(-Y))
          %>% mutate(fY = dbeta(Y, 2, 1),
                     fU = -2*log(U)/U)
) %>% head

          Y         U        fY       fU
1 0.6922911 0.5004282 1.3845822 2.766795
2 0.9627615 0.3818370 1.9255229 5.042788
3 0.7748989 0.4607504 1.5497978 3.363639
4 0.7604383 0.4674615 1.5208767 3.253480
5 0.5913213 0.5535954 1.1826425 2.136294
6 0.3574831 0.6994345 0.7149663 1.022206

Plotting results

ggplot(data = sim_df, aes(x = Y)) + 
    geom_histogram(aes(y=after_stat(density)), 
                   fill = 'goldenrod', color ='black',
                   center = 0.02, binwidth = 0.04)+
    geom_line(aes(y = fY), linewidth = 1) + 
    xlim(c(0,1)) + 
    labs(x = 'y', y = expression(f[Y](y))) + 
    theme_classic(base_size = 14)  +


ggplot(data = sim_df, aes(x = U)) + 
    geom_histogram(aes(y=after_stat(density)), 
                   fill = 'cornflowerblue',color ='black',
                   center = exp(-1)+.02,  binwidth = 0.04)+
    geom_line(aes(y = fU), linewidth = 1) + 
    xlim(c(0,1)) + 
    labs(x = 'u', y = expression(f[U](u))) + 
    theme_classic(base_size = 14)

The Pareto example

The Pareto distribution is named for the Italian Vilfredo Pareto
Originally used to model income inequality: a small percent of people own the greatest amount of wealth.
\(X\sim Pareto(\lambda,m)\) with \(\lambda >0\) and \(m>0\) if:

\[ f_X(x) = \left\{\begin{array}{ll} \frac{\lambda m^\lambda}{x^{\lambda+1}} & x > m \\ 0 & otherwise\\ \end{array} \right. \]

Pareto from exponential

If \(Y\sim EXP(\lambda)\), show that with \(m>0\), \(U = me^{Y}\sim Pareto(\lambda,m)\).

Support: \(m < U < \infty\)

\[g^{-1}(u) = \ln\left(\frac{u}{m}\right) \Rightarrow \frac{d}{du}g^{-1}(u) = \frac{m}{u}\cdot \frac{1}{m} = \frac{1}{u} \]

\[ \Rightarrow f_U(u) = f_Y(g^{-1}(u)) \cdot \frac{d}{du}g^{-1}(u) = \lambda e^{-\lambda \ln\left(\frac{u}{m}\right)}\cdot \frac{1}{u}\]

\[ = \lambda e^{\ln\left(\frac{u}{m}\right)^{-\lambda}}\cdot \frac{1}{u} = \lambda \frac{m^\lambda}{u^\lambda}\frac{1}{u} = \frac{\lambda m^\lambda}{u^{\lambda + 1}}, u > m\]

Verifying via simulation study

Verifying this result for a grid of \(\{\lambda, m\}\) values:

library(purrrfect)
N <- 10000

#Function to evaluate analytic pareto pdf:

dpareto <- \(u,lambda,m) {
  ifelse(u > m,(lambda*m^lambda)/(u^(lambda+1)),0)
}

# Simulating over a grid:

(pareto_sim <- parameters(~lambda, ~m,
           seq(2, 6, by = 2), c(1,2,4)
           )  
  %>% add_trials(N)
  %>% mutate(Y = map_dbl(lambda, .f = \(l) rexp(1, rate =l)))
  %>% mutate(U = m*exp(Y))
  %>% mutate(fU = dpareto(U, lambda, m))
) %>% head

# A tibble: 6 × 6
  lambda     m .trial      Y     U     fU
   <dbl> <dbl>  <dbl>  <dbl> <dbl>  <dbl>
1      2     1      1 1.66    5.27 0.0137
2      2     1      2 0.0593  1.06 1.67  
3      2     1      3 0.305   1.36 0.800 
4      2     1      4 0.655   1.92 0.280 
5      2     1      5 0.480   1.62 0.474 
6      2     1      6 0.461   1.59 0.502

Plotting the results

ggplot(data = pareto_sim) + 
    geom_histogram(aes(x = U,y=after_stat(density)), 
                   fill = 'goldenrod')+
    geom_line(aes(x = U, y = fU), col='cornflowerblue', size = 1) + 
    theme_classic(base_size = 16) + 
    scale_x_continuous(breaks = seq(0,14, by = 2), limits = c(0,15)) + 
    facet_grid(lambda ~ m, scales = 'free_y', labeller = label_both)

Simulated densities U with Pareto density superimposed