1
$\begingroup$

When trying to code this in R, I'm getting very confused about what to do. Apologies if my terminology is incorrect but I would be grateful for any advice.

The Problem:

I have been given two normal distributions for the mean and standard deviation of rat weights. So:

  1. The mean rate weight is itself a normal distribution with a mean of 1.68 and a standard deviation, or confidence interval, of 1.81 (mean = 12.68, SD = 1.81)
  2. The standard deviation is itself a normal distribution with a mean of 11.19 and a standard deviation of 3.2 (mean = 11.19, SD = 3.2)

Summary Mean Rate Weight: dnorm(Mean = 12.68, SD = 1.81) Standard Distribution of Rat Weight: dnorm(Mean = 11.19, SD = 3.2)

The Question:

In R, how to a code this to have a Monte Carlo run of 50,000 samples? Is the following example correct?

 MC_Runs = 50000

 Rat_Weight = rnorm(MC_Runs,
                   mean = rnorm(MC_Runs,mean = 12.68, sd = 1.81),
                   sd= rnorm(MC_Runs,mean = 11.19, sd = 3.2))
|cite|improve this question
$\endgroup$
  • 2
    $\begingroup$ This is a bit unusual as a standard deviation generally wouldn't be normal distributed since it is always positiv. $\endgroup$ – Roland Aug 27 '14 at 13:20
  • 2
    $\begingroup$ You should expect twelve failures of this code (on the average) because the generated SD will be negative: pnorm(-11.19/3.2) * 50000 is about $11.8$. Although that is not many, it should be taken seriously as a warning that this model might not adequately describe the rate weights and that the Monte Carlo output could have some misleading results in it. $\endgroup$ – whuber♦ Aug 27 '14 at 19:36
1
$\begingroup$

Your solution is correct, assuming the two normal random variables are independent. According to the R documentation of rnorm, you can input a vector of means and standard deviations for the mean and sd arguments respectively.

To verify, consider this toy example:

n <- 3 
mean_vector <- c(0,10,100)
sd_vector <- c(1,1,1)

rnorm(3, mean=mean_vector, sd=sd_vector)

Some output:

[1]  1.049676 11.566033 98.481899
[1] -1.374753  9.078215 99.465803
[1]  3.056377  9.837055 98.842553

Clearly the first variate for each simulation is $N(0,1)$ distributed, the second is $N(10,1)$ distributed, and the third is $N(100,1)$ distributed.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thank you @ Comp_Warrior and @Statsspecialist. I have been trying to verify this from first principal derivation but have had no luck. I wonder if you might know some source that I can look up for a reference? $\endgroup$ – hoof_hearted Aug 27 '14 at 12:56
  • $\begingroup$ @hoof_hearted I am unsure about what it is you are trying to verify - can you tell us? $\endgroup$ – Comp_Warrior Aug 27 '14 at 12:59
  • $\begingroup$ I was wondering if there is a mathematical derivation that will show me how the Mean Rate Weight [dnorm(Mean = 12.68, SD = 1.81)] and Standard Distribution of Rat Weight [dnorm(Mean = 11.19, SD = 3.2)] combine to give a final distribution. I want to be able to back up my code with a first principal check, but don't really know where to start. Hope this is clearer. $\endgroup$ – hoof_hearted Aug 27 '14 at 13:13
  • 1
    $\begingroup$ @hoof_hearted I think you are dealing with a conditional distribution here. So if $Z$ is the weight of the rats, $X$ is the distribution of the mean, and $Y$ is the distribution of the standard deviations, $Z|(X,Y)\sim N(X,Y)$. $X$ and $Y$ are independent normally distributed here. The joint distribution of the three would be $f_{X,Y,Z}(x,y,z) = f_{Z|(X,Y)}(z|x,y) f_{X,Y}(x,y)$. To find the (unconditional) distribution of $Z$, you would have to integrate the joint distribution over all possible values of $X$ and $Y$. $\endgroup$ – Comp_Warrior Aug 27 '14 at 13:23
  • 1
    $\begingroup$ Many thanks. I'll start down the road of conditional distributions and see where it gets me. Your help is much appreciated! $\endgroup$ – hoof_hearted Aug 27 '14 at 13:28
1
$\begingroup$

The code will work, and will give you random samples from a normal distribution with the set of means and SDs with the defined distributions. What you should note is that the generated sample of 50000 might not explore all the candidate pairs of the mean and SD values that are used to generate it. If that is of concern to you, then you might want to try Latin hyper-cube sampling.

|cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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