Sampling from a mixture of distributions (II)
The complete algorithm for sampling from a mixture distribution is:
- Choose a component.
- Generate a normal random number using the mean and standard deviation of the selected component.
choose_component()
, from the last exercise, is provided. Here you'll complete the second step and complete the definition of rmix()
.
This exercise is part of the course
Optimizing R Code with Rcpp
Exercise instructions
- Check that there are as many standard deviations as weights. That is, the size of
sds
is the same asd
. - Calculate
total_weight
as the sum of the weights. - Choose a component by calling
choose_component()
. - Simulate from the chosen component by generating a normal random number with the
j
th element ofmeans
andsds
.
Hands-on interactive exercise
Have a go at this exercise by completing this sample code.
#include
using namespace Rcpp;
// From previous exercise; do not modify
// [[Rcpp::export]]
int choose_component(NumericVector weights, double total_weight) {
double x = R::runif(0, total_weight);
int j = 0;
while(x >= weights[j]) {
x -= weights[j];
j++;
}
return j;
}
// [[Rcpp::export]]
NumericVector rmix(int n, NumericVector weights, NumericVector means, NumericVector sds) {
// Check that weights and means have the same size
int d = weights.size();
if(means.size() != d) {
stop("means size != weights size");
}
// Do the same for the weights and std devs
if(___) {
stop("sds size != weights size");
}
// Calculate the total weight
double total_weight = ___;
// Create the output vector
NumericVector res(n);
// Fill the vector
for(int i = 0; i < n; i++) {
// Choose a component
int j = ___(___, ___);
// Simulate from the chosen component
res[i] = ___::___(___, ___);
}
return res;
}
/*** R
weights <- c(0.3, 0.7)
means <- c(2, 4)
sds <- c(2, 4)
rmix(10, weights, means, sds)
*/