0

I am ask to sample and then do a marginal distribution $f(x_1,x_5)$ and plot it. I have the following code which works but dnorm is used for one dimension so I was wondering if i need to change it to dmvnorm.

If so, i change mu=mu.marginal, sigma=sigma.marginal, added a y sample, but dmvnorm says error because of non array input. Anye help is appreciated.

Model of multivariable normal:

mu = c(1,2,6,2,-4)
sigma = c(1,2,1,0.5,2)
rho = diag(rep(1,5))
rho[1,2] = rho[2,1] = 0.4
rho[1,3] = rho[3,1] = -0.3
rho[1,4] = rho[4,1] = -0.7
rho[3,5] = rho[5,3] = 0.2
rho[4,5] = rho[5,4] = 0.5
Sigma = rho * (sigma %o% sigma)

my code:

p = c(1,5)
(mu.marginal = mu[p])
(Sigma.marginal = Sigma[p,p])
# p is one-dimensional: use dnorm() to compute marginal pdf
x = seq(-1,6,by=0.01)

fx = dnorm(x,mean=mu.marginal,sd=sqrt(Sigma.marginal))
ggplot(data=data.frame(x=x,y=fx),mapping=aes(x=x,y=y)) + geom_line(col="blue")
Isabelle
  • 7
  • 2
  • Your code comment says "p is one-dimensional", but `(x_1, x_5)` is two-dimensional... To make sure I understand properly your task, you are to sample from the marginal distribution of `(x_1, x_5)` and plot the samples? – duckmayr Jul 24 '20 at 03:14
  • yes sorry its meant to be p is two dimensional – Isabelle Jul 24 '20 at 04:20

1 Answers1

0

It seems to me you were on the right track with mvtnorm and came close to a solution... I'm not sure how you ran into a non-array input error, but here's what I got with using mvtnorm:

set.seed(123)
dat <- mvtnorm::rmvnorm(1e4, mean = mu.marginal, sigma = Sigma.marginal)
dat <- as.data.frame(dat)

ggplot(dat, aes(x = V1, y = V2)) +
    geom_bin2d()

enter image description here

You can see it's fairly spherical, which is what you'd expect since the off-diagonal elements of Sigma.marginal are 0 (which means that x_1 and x_5 are marginally independently normally distributed...)

duckmayr
  • 16,303
  • 3
  • 35
  • 53
  • thank you, what would cause array errors usually in this case? ill try to replicate the code that yielded my error if i can – Isabelle Jul 24 '20 at 04:20
  • i used dmvnorm instead of mvtnorm::rmvnorm is there a reason why it might have happend? – Isabelle Jul 24 '20 at 04:21
  • @Isabelle Yes, perhaps; it depends on how specifically you called it (i.e. what arguments you gave. Without seeing the code you used to call `dmvnorm()`, my best guess would be that you passed something that was not a matrix as `sigma` – duckmayr Jul 24 '20 at 10:41