Final
Moiya Josephs
2023-04-20
# Load the required library for Gamma distribution
library(stats)Problem 1
Problem Density 1: X~Gamma
Using R, generate a random variable X that has 10,000 random Gamma pdf values. A Gamma pdf is completely described by n (a size parameter) and lambda \(\lambda\) (a shape parameter). Choose any n greater 3 and an expected value \(\lambda\) between 2 and 10.
set.seed(0)
# Parameters for Gamma distribution
n <- 5 # size parameter
lambda <- 5 # shape parameter
# Generate random values from Gamma distribution
X <- rgamma(10000, shape = n, rate = lambda)
# Plot the histogram of the generated data
hist(X, freq = FALSE, col = "lightblue", main = "Gamma Distribution", xlab = "X", ylab = "Density")
# Overlay the theoretical Gamma pdf on the histogram
curve(dgamma(x, shape = n, rate = lambda), add = TRUE, col = "red", lwd = 2)
The rgamma function from the stats library is used to
generate 10,000 random values from a Gamma distribution with the number
of values to generate set to n and shape set to lambda.
Problem Density 2: Sum of Exponentials
Here I generate and plot random variable Y, which is the sum of Exponentials. It should also be noted that the sum of exponentials is also a gamma distribution.The histrogram below is very similar to the X histogram above. That is because the shape and size parameter are the same.
set.seed(0)
Y<- rexp(10000,5)+rexp(10000,5)+rexp(10000,5)+rexp(10000,5)+rexp(10000,5)
hist(Y, freq = FALSE, col = "lightblue", main = "Sum of Exponentials Distribution", xlab = "Y", ylab = "Density")
curve(dgamma(x, shape = n, rate = lambda), add = TRUE, col = "red", lwd = 2)
Probability Density 3: Exponential
For Probability density function is a single exponential with the same lambda as the previous pdfs.
# set the rate/shape parameter
set.seed(0)
# generate 10,000 observations from the exponential distribution
Z <- rexp(n=10000, rate=lambda)
hist(Z, freq = FALSE, col = "lightblue", main = "Exponential Distribution", xlab = "Z", ylab = "Density")
curve(dexp(x, rate = lambda), add = TRUE, col = "red", lwd = 2)
1b
Calculate the empirical expected value (means) and variances of all three pdfs.
X Gamma
The mean of gamma distribution is given by the product of its shape and scale parameters:
The formula for the mean of a gamma distribution is:
\[ \frac{\alpha}{\beta} \] Where \(\alpha\) is the shape parameter and \(\beta\) is the rate.
For X both parameters are 5, so the mean should be equal to 1. This
can be tested as shown below, by applying the formula or calling the R
function mean on the random variable X.
# shape a, rate b
# n/ lambda
n/lambda## [1] 1mean(X)## [1] 1.000132Variance of a gamma distribution is:
\[ Var(X) = \frac{\alpha}{\beta^2} \]
n/lambda^2## [1] 0.2var(X)## [1] 0.1993912Y Sum of Exponentials
Theoretical Mean
As discussed above Y is a gamma distribution. The mean of gamma distribution is given by the product of its shape and scale parameters:
Where shape is 5 and the scale is also 5.
\[ \frac{1 * 5}{5} = 1 \]
Empirical Mean
mean(Y)## [1] 1.002137Variance
The variance of a gamma distribution is:
\[ Var(Y) = \frac{\alpha}{\beta^2} \]
5/5^2## [1] 0.2var(Y)## [1] 0.2070123Z Exponential
Emperical Mean for an Exponential
\[ E(X)= 1/ \lambda \]
So for this we are expecting \(1/5\).
mean(Z)## [1] 0.20121781/5## [1] 0.2Variance
\[ E(X)= 1/ \lambda^2 \]
1/5^2## [1] 0.04var(Z)## [1] 0.04194312Problem 1b:
Using calculus, calculate the expected value and variance of the Gamma pdf(X). Using the moment generating function for exponentials, calculate the expected value of a single exponential (Z) and the sum of the exponentials.
Expected value and variance of X (Using Calculus)
The probability density function of a Gamma distribution is given below: (Source: https://proofwiki.org/wiki/Expectation_of_Gamma_Distribution)
\[ f_x(x)=\frac{\beta^\alpha x^{\alpha-1}e^{-\beta x}}{\gamma(\alpha)} \] The expected value is the integral with the bounds from 0 to infinity of the function \(x*f_x(x)\). We could do it mathematically to get the result:
\[ E(X) = \int_0^\infty xf_x(x)dx=\frac{\alpha}{\beta} \]
Where in this problem \(\alpha = 5\) and \(\beta = 5\).
We can also use R to find the integral of the function like in this Source.
integral.X <- function(x) x * dgamma(x, n, lambda)
exp.X <- integrate(integral.X, lower = 0, upper = Inf)$valueThe variance is given by the following: (Source: https://proofwiki.org/wiki/Variance_of_Gamma_Distribution)
\[ Var(X) = \int_0^\infty x^2 f_x(x)dx - (E(X))^2=\frac{\alpha}{\beta^2} \]
In the context of this problem, \(\alpha = 5\) and \(\beta = 5\) so the variance should be equal to \(\frac{5}{25}\).
5/25## [1] 0.2integral.X2 <- function(x) x^2 * dgamma(x, n, lambda)
varX <- integrate(integral.X2, lower = 0, upper = Inf)$value
varX - exp.X^2## [1] 0.2Expected value and variance of Z (Using Moment Generating function)
The moment generating function for an exponential distribution is defined below. This is also the expected value. Thorugh calculus it is deduced to be \(\frac{1}{\lambda}\).
\[ M_X(t) = E[e^{tX}]\\ = \int_{0}^{\infty}xp_X(x)dx \\ = \lambda\int_{0}^{\infty}e^{tx}e^{-x(\lambda)}dx\\ = \lambda [\frac{1}{(t-\lambda)}e^{x(t-\lambda)}]_{0}^{\infty}\\ = \frac{\lambda}{\lambda-t}\\ \text{When } t=0, \frac{1}{\lambda} \] Here we evaluate for when t = 0
1/lambda## [1] 0.2As E(X) is \(1/lambda\) the variance is \(E(X^2)-(E(X))^2\).
OR using R to integrate:
integral.Z <- function(x) x * dexp(x, n)
exp.Z <- integrate(integral.Z, lower = 0, upper = Inf)$value
exp.Z## [1] 0.2The variance:
\[ E(X^2)=M_X^{(2)}(t)= \frac{(2\lambda)}{\lambda - t}|_{t=0} = \frac{2}{\lambda^2} \] \[ E(X^2)-(E(X))^2 \\ = \frac{2}{\lambda^2} - (\frac{1}{\lambda})^2\\ =\frac{1}{\lambda^2} \] Which is the same as above.
1/lambda^2## [1] 0.04Expected value and variance of Y (Using Moment Generating function)
Y is the sum of same exponential distributions n times. So we can multiply the moment generating function by n times to get the mean and variance.
\[ 5 * M_Z(0)\\ =5*\frac{1}{5} \]
The variance will be the same as defined before:
\[ 5 * M^2_Z(0) = 5 * (1/5^2) = 5/25 \]
Problem 1 c-e
For PDF Z, calculate emperically probabilities a through c. Then evaluate through calculus if the memoryless property holds.
The PDF for the exponential distribution is defined below.
\[ f(x)= \frac{1}{\mu}e^{-\frac{1}{\mu}}x \]
To calculate probabilities for specific PDF, we have to use the CDF. The CDF is the integral of the PDF.
\[ F(x)=\int_0^\infty[\frac{1}{\mu}e^{-\frac{1}{\mu}}x]=1-e^{-\frac{x}{\mu}} \]
A
\[ P(Z > \lambda | Z> \lambda/2) \]
\[ \lambda =5 \]
\[ P(Z >5 | Z>5/2) = P(Z>5|Z>2.5) \]
\[ P(Z>5|Z>2.5)\\= \frac{1-P(Z<5)}{1-P(Z<2.5)}\\= \frac{1-F(5)}{1-F(2.5)}\\=\frac{1-(1-e^{-5/5})}{1-(1-e^{-2.5/5})}\\=\frac{e^{-1}}{e^{-1/2}}\\=e^{-1-(-1/2)}\\=e^{-1/2} \]
Given the memoryless property:
\[ P(Z > r + t| Z>r) = P(Z>t)\\ r = 2.5, t=2.5 \]
\[ P(Z > 2.5 + 2.5| Z> 2.5) = P(Z>2.5)\\ P(Z>2.5)=1-P(Z<2.5)=1-(1-e^{-1/5*(2.5)})\\=1-1+e^{-0.5} = \\ e^{-0.5} \]
exp(-.5)## [1] 0.6065307Another way to test this answer is to use the pexp
function in R. Where the pexp gives us the \(P(Z<\lambda)\) so to find \(P(Z<\lambda)\) we take the complement
like we did above.
(Z.a=1-pexp(2.5,(1/5)))## [1] 0.6065307The emperical probability is: 0.6065307
The following calculus proof for the above probability.
The exponential distribution probability density function is:
\[ f(x) = \lambda e^{-\lambda x} \]
To calculate the empirical probability \(P(Z > \lambda | Z > \lambda/2)\), we need to use the definition of conditional probability:
\[ P(Z > \lambda | Z > \lambda/2) = \frac{P(Z > \lambda \cap Z > \lambda/2)}{P(Z > \lambda/2)} \]
Also useful is the rule for intersection:
\[ P(A\cap B)={P(A/B)}{P(B)} \]
As for the memoryless property it states,
\[ P(Z > r + t| Z>r) = P(Z>t)\\ \]
In context with this problem, \(r+t = \lambda\) and \(r=\lambda/2\). Because \(lambda=5\) and \(r=2.5\) we know in order for \(r+t=5\), \(t=2.5\).
\[ =\frac{P(Z > \lambda \cap Z > \lambda/2)}{P(Z > \lambda/2)}\\ P(Z > 2.5 + 2.5 | Z > 2.5) = \\ \frac{P(Z > 5 \cap Z > 2.5)}{P(Z > 2.5)} \\ =P(Z>2.5) \]
We can use the cumulative distribution function (cdf) to calculate the probabilities:
\[ \int_\lambda^\infty e^{-\lambda x}dx \]
To calculate empirically using this formula we get:
For \(P(Z>5)\)
\[ \int_5^\infty e^{5x}dx \\ =e^{-5x} \]
pdf <- function(x, lambda) exp(-x/lambda)/lambda(Z.p1 <- integrate(pdf, 5, Inf, lambda)$value)## [1] 0.3678794For \(P(Z>2.5)\)
\[ \int_{2.5}^\infty e^{-2.5x}dx \\ =e^{-2.5x} \]
which in R would be:
(Z.p2 <- integrate(pdf, 2.5, Inf, lambda)$value)## [1] 0.6065307Now we need to apply the intersection rule along with the conditional probability rule. If the function is memoryless, we should get the same as \(P(Z>2.5)\)
((Z.p1*Z.p2)/Z.p2)/Z.p2## [1] 0.6065307We get the same result so it proves the memoryless property holds.
–
B
\[ P(Z >2\lambda | Z>\lambda) = P(Z>5(2)|Z>5)\\=P(Z>10|Z>5) \]
\[ P(Z > 5 + 5| Z> 5) = P(Z>5)\\ P(Z>5)=1-P(Z<5)=1-(1-e^{-1/5*5})\\=1-1+e^{-1} \]
exp(-1)## [1] 0.3678794(Z.b=1-pexp(5,(1/5)))## [1] 0.3678794The emperical probability is:
\[ P(Z >2\lambda | Z>\lambda) = P(Z>5(2)|Z>5)\\=P(Z>10|Z>5) = P(Z>5) \]
For \(P(Z>10)\)
\[ \int_{10}^\infty e^{-10x}dx \\ =e^{-10x} \]
(Z.p1 <- integrate(pdf, 10, Inf, lambda)$value)## [1] 0.1353353For \(P(Z>5)\) we can use R again
to find the value and save it into Z.p2.
(Z.p2 <- integrate(pdf, 5, Inf, lambda)$value)## [1] 0.3678794We should get the same result as Z.p2.
((Z.p1*Z.p2)/Z.p2)/Z.p2## [1] 0.3678794We get the same result so it proves the memoryless property holds.
C
\[ P(Z >3\lambda | Z>\lambda) = P(Z>5(3)|Z>5)\\=P(Z>15|Z>5) \]
This follows that:
\[ r+t = 15\\ r = 5 \\ t = 10 \]
So we get:
\[ P(Z > 5 + 10| Z> 5) = P(Z>10)\\ P(Z>10)=1-P(Z<10)=1-(1-e^{-1/5*10})\\=1-1+e^{-2} \]
exp(-2)## [1] 0.1353353(Z.c=1-pexp(10,(1/5)))## [1] 0.1353353The empirical probability is:
\[ P(Z >3\lambda | Z>\lambda) = P(Z>5(3)|Z>5)\\=P(Z>15|Z>5) = P(Z>10) \]
(Z.p1 <- integrate(pdf, 15, Inf, lambda)$value)## [1] 0.04978707(Z.p2 <- integrate(pdf, 5, Inf, lambda)$value)## [1] 0.3678794(Z.p3 <- integrate(pdf, 10, Inf, lambda)$value)## [1] 0.1353353((Z.p1*Z.p2)/Z.p2)/Z.p2## [1] 0.13533531d. Marginal and joint probabilities
\[ P(X \cap Y) = P(X|Y) \times P(Y) \]
To make a marginal probability, first you multiply the random values together. Then you sum the columns and then the rows to get the total marginal distribution. To find the joint probability the total sum is used as a divisor and the individual product from YZ, the dividend.
I perform these steps below.
The quantiles for the Y or Z distribution are defined below and put
Y.quant and Z.quant.
(Y.quant <- quantile(Y, probs=c(.25,.50,.75,1)))## 25% 50% 75% 100%
## 0.6698686 0.9367503 1.2606274 3.4414941(Z.quant <- quantile(Z, probs=c(.25,.50,.75,1)))## 25% 50% 75% 100%
## 0.05704276 0.13958030 0.27605182 1.83689029I converted the variables to matrixes for easy multiplication
(Z.mat <- matrix(Z.quant,nrow=4,ncol=1,byrow=TRUE))## [,1]
## [1,] 0.05704276
## [2,] 0.13958030
## [3,] 0.27605182
## [4,] 1.83689029(Y.mat <- matrix(Y.quant,nrow=1,ncol=4,byrow=TRUE))## [,1] [,2] [,3] [,4]
## [1,] 0.6698686 0.9367503 1.260627 3.441494The product is them stored into ZY and the
(ZY <- Z.mat %*% Y.mat)## [,1] [,2] [,3] [,4]
## [1,] 0.03821115 0.05343482 0.07190967 0.1963123
## [2,] 0.09350046 0.13075188 0.17595875 0.4803648
## [3,] 0.18491843 0.25859161 0.34799848 0.9500307
## [4,] 1.23047506 1.72070746 2.31563423 6.3216471By summing all columns and rows we get the below two way frequency table.
ZY_marginal <- rbind(ZY,colSums(ZY))
ZY_marginal <- cbind(ZY_marginal,rowSums(ZY_marginal))
ZY_marginal## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.03821115 0.05343482 0.07190967 0.1963123 0.3598680
## [2,] 0.09350046 0.13075188 0.17595875 0.4803648 0.8805759
## [3,] 0.18491843 0.25859161 0.34799848 0.9500307 1.7415392
## [4,] 1.23047506 1.72070746 2.31563423 6.3216471 11.5884639
## [5,] 1.54710511 2.16348577 2.91150113 7.9483550 14.5704470The joint probability is done below, where I divide each cell by the total sum of the matrix.
ZY_joint<-t(apply(ZY,1, function(x) x/sum(ZY)) )
ZY_joint## [,1] [,2] [,3] [,4]
## [1,] 0.002622511 0.003667343 0.00493531 0.01347332
## [2,] 0.006417130 0.008973773 0.01207641 0.03296843
## [3,] 0.012691336 0.017747679 0.02388386 0.06520258
## [4,] 0.084450056 0.118095722 0.15892678 0.43386776ZY_joint <- rbind(ZY_joint,colSums(ZY_joint))
(ZY_joint <- cbind(ZY_joint,rowSums(ZY_joint)))## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.002622511 0.003667343 0.00493531 0.01347332 0.02469849
## [2,] 0.006417130 0.008973773 0.01207641 0.03296843 0.06043575
## [3,] 0.012691336 0.017747679 0.02388386 0.06520258 0.11952545
## [4,] 0.084450056 0.118095722 0.15892678 0.43386776 0.79534032
## [5,] 0.106181033 0.148484516 0.19982236 0.54551209 1.00000000The final joint probability shows distribution for each Y and Z combination and the total probability equals 1. If the probabilities are independent then the conditional probability can be determined by Y and Z.
To see if this can explain P(YZ)=P(Y)P(Z). If it is
independent we can multiply the probabilities of each random variable to
get the conditional probability.
We can check for independence by:
\[ f(x,y)=f_X(x)f_Y(y) \]
Where \(f\) is the contingency table, or the marginal table.
\[ P(Z=1,Y=1) = P(Z=1)*P(Y=1)\\ ZY[1,1] = sum(Z[1,])sum(Y[,1])\\ \]
We will choose Z[1,] and Y[,1] which is
ZY[1,1]## [1] 0.03821115This should be equal to the sum of the Z[1,] and Y[,1] which is given in the marginal table above. By using R we get.
ZY_marginal[1,5] * ZY_marginal[5,1]## [1] 0.5567536As we can see the two are not equal.
1e. Check to see if independence holds
Z.Y<-(Z.mat%*%Y.mat)Z.Y## [,1] [,2] [,3] [,4]
## [1,] 0.03821115 0.05343482 0.07190967 0.1963123
## [2,] 0.09350046 0.13075188 0.17595875 0.4803648
## [3,] 0.18491843 0.25859161 0.34799848 0.9500307
## [4,] 1.23047506 1.72070746 2.31563423 6.3216471Z.Y<- rbind(Z.Y,colSums(Z.Y))
Z.Y## [,1] [,2] [,3] [,4]
## [1,] 0.03821115 0.05343482 0.07190967 0.1963123
## [2,] 0.09350046 0.13075188 0.17595875 0.4803648
## [3,] 0.18491843 0.25859161 0.34799848 0.9500307
## [4,] 1.23047506 1.72070746 2.31563423 6.3216471
## [5,] 1.54710511 2.16348577 2.91150113 7.9483550Z.Y<- cbind(Z.Y,rowSums(Z.Y))
Z.Y## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.03821115 0.05343482 0.07190967 0.1963123 0.3598680
## [2,] 0.09350046 0.13075188 0.17595875 0.4803648 0.8805759
## [3,] 0.18491843 0.25859161 0.34799848 0.9500307 1.7415392
## [4,] 1.23047506 1.72070746 2.31563423 6.3216471 11.5884639
## [5,] 1.54710511 2.16348577 2.91150113 7.9483550 14.5704470Fisher’s Exact Test and Chi Square Test are both statistical tests used to determine the independence between two categorical variables in a contingency table.
The hypotheses we are testing are:
\[ H_0: \text{the variables are independent} \\ H_1: \text{the variables are dependent} \]
Chi Square
The Chi Square Test assumes that the expected cell frequencies are large enough for the distribution to be approximated.
(test<-chisq.test(Z.Y))## Warning in chisq.test(Z.Y): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: Z.Y
## X-squared = 3.9367e-31, df = 16, p-value = 1This gives us a warning because the smallest expected frequencies are lower than 5. To avoid this error we will use the Fisher’s exact test which does not have the assumption of the minimum expected value.
Fisher’s Exact Test a test that calculates the exact probability of obtaining the observed frequencies in a contingency table, given the marginals (row and column totals). It is useful when dealing with small sample sizes or when expected cell frequencies are less than 5.
(test2<-fisher.test(Z.Y))## Warning in fisher.test(Z.Y): 'x' has been rounded to integer: Mean relative
## difference: 0.095473##
## Fisher's Exact Test for Count Data
##
## data: Z.Y
## p-value = 1
## alternative hypothesis: two.sidedtest2$p.value## [1] 1So we fail to reject the null hypothesis.
The most appropriate test depends on the ssize of the data. Fisher’s Exact Test is more appropriate when the sample size is small and the expected cell frequencies are less than 5, or when the assumptions of the Chi Square Test are not met. On the other hand, the Chi Square Test can be used when dealing with larger sample sizes and expected cell frequencies are greater than 5.