Can any Probability Distribution be Integrated?

Question

Can any Probability Distribution be Integrated?

Suppose I have a Mixture Model. Mixture Models can be considered as a "linear combination of probability distribution functions". Combining several probability distribution functions together allows Mixture Models to "explain" complex patterns within the data that a single probability distribution would not be able to explain as easily. For example, different Normal Distributions can be combined together to make a Mixture Model (note: contrary to common logic, a mixture of two Normal Distributions will not necessarily result in a Normal Distribution):

Suppose I want to create the following Mixture Model by combining two Normal Distributions:

**P(x,y,z) = 0.3 * N1 + 0.7 * N2

N1(x,y,z) ~ Normal with mu = (1,1,1) and sigma(1,0,0,0,1,0,0,0,1)

N2(x,y,z) ~ Normal with mu = (1,0,1) and sigma(1,0,0,0,1,0,0,0,1)**

Question: If I fix the value of x = 2, can I use Monte Carlo Integration to integrate this conditional probability distribution? (over some arbitrary range)

I tried to demonstrate this using the R programming language.

Part 1: I defined the conditional distribution of this mixture model for x = 2:

#define constants needed for the multivariate normal
    sigma1 <- c(1,0,0,0,1,0,0,0,1)
      sigma <- matrix(sigma1, nrow=3, ncol= 3, byrow = TRUE)
      sigma_inv <- solve(sigma)
      sigma_det <- det(sigma)
      denom = sqrt( (2*pi)^4 * sigma_det) 

#mixture model
    target <- function(y,z)
      
    {
      x_one = 2 - 1
      x_two = y - 1
      x_three = z - 1
     
        
      x_t = c(x_one, x_two, x_three)
      x_t_one <- matrix(x_t, nrow=3, ncol= 1, byrow = TRUE)
      x_t_two =  matrix(x_t, nrow=1, ncol= 3, byrow = TRUE)
      
      
     
      num = exp(-0.5 * x_t_two  %*%  sigma_inv  %*%  x_t_one)
        
      answer_1 = num/denom
 
    
      x_one2 = 2 - 1
      x_two2 = y - 1
      x_three2 = z - 1
     
        
      x_t2 = c(x_one2, x_two2, x_three2)
      x_t_one2 <- matrix(x_t2, nrow=3, ncol= 1, byrow = TRUE)
      x_t_two2 =  matrix(x_t2, nrow=1, ncol= 3, byrow = TRUE)
      
      
      # In this part, as it's (x-mu)^T * SIGMA * (x-mu)
      
   
      num2 = exp(-0.5 * x_t_two2  %*%  sigma_inv  %*%  x_t_one2)
        
      answer_2 = num2/denom
      return(0.3*answer_1 + 0.7*answer_2)
    }

Part 2: I then attempted to integrate this conditional probability distribution using Monte Carlo Integration:

means_vec<-c()
for(j in 1:10000){
    set.seed(j)
    u_y<-runif(10000)
    u_z<-runif(10000)
    sim_vec<-c()
    for(i in 1:length(u_y)){
        sim_vec[i]<-target(y=u_y[i],z=u_z[i])
    }
    means_vec[j]<-mean(sim_vec)
}
print(means_vec)

The above code took 10,000 random points (from a uniform distribution) for "y" and "z", and evaluated the function (i.e. the conditional mixture distribution) at these 10,000 random points : then, using the Monte Carlo Estimator, the integral of this conditional normal distribution was calculated (by averaging the values of the function evaluated at these 10,000 random points). Then, this entire process is repeated 100 times : the integral of the conditional distribution can be considered as the average integral from each of these 100 iterations:

head(means_vec)

[1] 0.01124573 0.01123975 0.01123244 0.01123575 0.01124451 

#final answer of the integral
mean(means_vec)

[1] 0.01125189

Can someone please tell me if what I have done is correct? In practice, is this how the conditional distribution of a mixture models is integrated?

Answer 1

There are four little errors.

• In computing the denominator you take a power of 4 (instead of 3)

• You had a wrong mean vector x_two2 = y - 1

• You integrate over the joint distribution f(x,y,z) instead of the conditional distribution f(x,y,z)/f(x)

  In the corrected code below it is simply dividing by dnorm(1) which can be done because x is independent from y and z.

• The range of integration is small (the runif functions)

#define constants needed for the multivariate normal
    sigma1 <- c(1,0,0,0,1,0,0,0,1)
      sigma <- matrix(sigma1, nrow=3, ncol= 3, byrow = TRUE)
      sigma_inv <- solve(sigma)
      sigma_det <- det(sigma)
      denom = sqrt( (2*pi)^3 * sigma_det)

#mixture model
    target <- function(y,z)
      
    {
      x_one = 2 - 1
      x_two = y - 1
      x_three = z - 1
     
        
      x_t = c(x_one, x_two, x_three)
      x_t_one <- matrix(x_t, nrow=3, ncol= 1, byrow = TRUE)
      x_t_two =  matrix(x_t, nrow=1, ncol= 3, byrow = TRUE)
      
      
     
      num = exp(-0.5 * x_t_two  %*%  sigma_inv  %*%  x_t_one)
        
      answer_1 = num/denom
 
    
      x_one2 = 2 - 1
      x_two2 = y - 0
      x_three2 = z - 1
     
        
      x_t2 = c(x_one2, x_two2, x_three2)
      x_t_one2 <- matrix(x_t2, nrow=3, ncol= 1, byrow = TRUE)
      x_t_two2 =  matrix(x_t2, nrow=1, ncol= 3, byrow = TRUE)
      
      
      # In this part, as it's (x-mu)^T * SIGMA * (x-mu)
      
   
      num2 = exp(-0.5 * x_t_two2  %*%  sigma_inv  %*%  x_t_one2)
        
      answer_2 = num2/denom
      return(0.3*answer_1 + 0.7*answer_2)
    }

means_vec<-c()
for(j in 1:10){
    set.seed(j)
    box_length = 10
    u_y<-runif(10000,-box_length/2,box_length/2)
    u_z<-runif(10000,-box_length/2,box_length/2)
    sim_vec<-c()
    for(i in 1:length(u_y)){
        sim_vec[i]<-target(y=u_y[i],z=u_z[i])
    }
    means_vec[j]<-mean(sim_vec)
}
print(means_vec*box_length^2/dnorm(1))

mean(means_vec)*box_length^2/dnorm(1)

Answer 2

All probability distributions have to integrate/sum to one by definition. If that integration fails, it is not a probability distribution.