2

I'm trying to create a Gaussian Mix function according to these parameters:

  • For each sample, roll a die with k sides
  • If the j-th side appears from the roll, draw a sample from Normal(muj, sdj) where muj and sdj are the mean and standard deviation for the j-th Normal distribution respectively. This means you should have k different Normal distributions to choose from. Note that muj is the mathematical form of referring to the j-th element in a vector called mus.
  • The resulting sample from this Normal is then from a Gaussian Mixture.

Where:

  • n, an integer that represents the number of independent samples you want from this random variable
  • mus, a numeric vector with length k
  • sds, a numeric vector with length k
  • prob, a numeric vector with length k that indicates the probability of choosing the different Gaussians. This should have a default to NULL.

This is what I came up with so far:

n <- c(1)
mus <- c()
sds <- c()
prob <- c()

rgaussmix <- function(n, mus, sds, prob = NULL){
  if(length(mus) != length(sds)){
    stop("mus and sds have different lengths")
  }
  for(i in 1:seq_len(n)){
    if(is.null(prob)){
      rolls <- c(NA, n)
      rolls <- sample(c(1:length(mus)), n, replace=TRUE)
      avg <- rnorm(length(rolls), mean=mus[rolls], sd=sds[rolls])
    }else{
      rolls <- c(NA, n)
      rolls <- sample(c(1:length(mus), n, replace=TRUE, p=prob))
      avg <- rnorm(length(rolls), mean=mus[rolls], sd=sds[rolls])
    }
  }
  return(avg)
}

rgaussmix(2, 1:3, 1:3)

It seems to match most of the requirements, but it keeps giving me the following error:

numerical expression has 2 elements: only the first usednumber of items to replace is not a multiple of replacement length

I've tried looking at the lengths of multiple variables, but I can't seem to figure out where the error is coming from!

Could someone please help me?

1 Answers1

0

If you do seq_len(2) it gives you:

[1] 1 2

And you cannot do 1:(1:2) .. it doesn't make sense

Also you can avoid the loops in your code, by sampling the number of tries you need, for example if you do:

rnorm(3,c(0,10,20),1)
[1] -0.507961  8.568335 20.279245

It gives you 1st sample from the 1st mean, 2nd sample from 2nd mean and so on. So you can simplify your function to:

rgaussmix <- function(n, mus, sds, prob = NULL){
  if(length(mus) != length(sds)){
    stop("mus and sds have different lengths")
  }
  if(is.null(prob)){
    prob = rep(1/length(mus),length(mus)) 
  }
  rolls <- sample(length(mus), n, replace=TRUE, p=prob)
  avg <- rnorm(n, mean=mus[rolls], sd=sds[rolls])
  avg
}

You can plot the results:

plot(density(rgaussmix(10000,c(0,5,10),c(1,1,1))),main="mixture of 0,5,10")

enter image description here

StupidWolf
  • 45,075
  • 17
  • 40
  • 72