0
$\begingroup$

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):

enter image description here

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?

$\endgroup$
10
  • 5
    $\begingroup$ All probability distributions have to integrate/sum to one over their support by definition. If that integration fails, it is not a probability distribution. $\endgroup$
    – Galen
    Nov 9, 2021 at 17:48
  • 1
    $\begingroup$ I don't have the time to check the details of how you're setting it up, but prima facie I'd say Monte Carlo integration is worth looking at if the symbolic form is difficult. $\endgroup$
    – Galen
    Nov 9, 2021 at 17:51
  • 1
    $\begingroup$ @Galen Not all functions with non-zero areas under their curves can be (analytically) integrated: this is one reason why numerical methods exist. $\endgroup$
    – Alexis
    Nov 9, 2021 at 19:18
  • 1
    $\begingroup$ @Alexis That's fine. Being able to write a closed form solution doesn't mean the integral doesn't exist. By definition the definite integral over the support of the function must exist, and be equal one, for it to be a distribution. $\endgroup$
    – Galen
    Nov 9, 2021 at 19:32
  • 1
    $\begingroup$ @Galen We agree. :) $\endgroup$
    – Alexis
    Nov 9, 2021 at 19:47

2 Answers 2

3
$\begingroup$

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)
$\endgroup$
11
  • $\begingroup$ Thank you so much for your answer! I really appreciate it! Can you please explain the significance of this line of code:.... print(means_vec*box_length^2/dnorm(1)) ...thank you so much! $\endgroup$
    – stats_noob
    Nov 10, 2021 at 6:26
  • 1
    $\begingroup$ @stats555 It contains two parts. (1) is a multiplication by the array over which you generate the random points. If you would do a grid integration you would multiply by the area divided by the number of points. For example n=100; x = seq(-4,4,8/n) ; sum(dnorm(x))*8/n; or like you do put the division by n in mean(dnorm(x))*8. The MC integration contains the same factor but the points are not in a regular grid. (2) The dnorm(1) is to make the function a conditional distribution. You need to compute $f(x,y,z|x) = f(x,y,z)/f(x)$ and before you only computed $f(x,y,z)$. $\endgroup$ Nov 10, 2021 at 6:43
  • $\begingroup$ @ Sextus Empiricus: Thank you for your reply! In my question, I was fixing the value of "x" at "x = 2". Therefore, should the last line of code be: mean(means_vec)*box_length^2/dnorm(2)? Thank you! $\endgroup$
    – stats_noob
    Nov 15, 2021 at 4:59
  • $\begingroup$ is this line of code optional? print(means_vec*box_length^2/dnorm(1)) .... that is, the only the final line of code is required? mean(means_vec)*box_length^2/dnorm(1) $\endgroup$
    – stats_noob
    Nov 15, 2021 at 5:00
  • $\begingroup$ If I were to calculate f(x,y,z|x,y) = f(x,y,z) / f(x,y) ... would this now require the density of a bivariate normal distribution in terms of "x and y"? Thank you so much for your help! $\endgroup$
    – stats_noob
    Nov 15, 2021 at 5:01
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.