CODEWITHSUNDEEP
r
adam.888
adam.888

Reputation: 7846

R obtaining a probability distribution

I have a relationship: y = a + b + c

I have the average and standard deviation of a, b and c and I would like to obtain the probability distribution of y from this by Monte Carlo simulation.

Is there a function or package or easy way that I can use to do this?

Upvotes: 1

Views: 294

Answers (2)

IRTFM
IRTFM

Reputation: 263332

There would be two approaches: bootstrapping which I think is what you might mean by MonteCarlo or if you are more interested in the theory than constructing estimates from empiric distributions, the 'distr' package and its friends 'distrSim" and "distrTEst".

require(boot)
ax <- rnorm(100); bx<-runif(100); cx<- rexp(100)
dat <- data.frame(ax=ax,bx=bx,cx=cx)

boot(dat, function(d){ with(d, mean(ax+bx+cx) )}, R=1000,  sim="parametric")
boot(dat, function(d){ with(d, sd(ax+bx+cx) )}, R=1000,  sim="parametric")

Upvotes: 2

Seth
Seth

Reputation: 4795

I assume that your are assuming your inputs a,b and c are normally distributed because you say you can define them with mean and standard deviation. If that is the case, you can do this pretty fast without any special package.

 mu.a=33
 mu.b=32
 mu.c=13
 sigma.a=22
 sigma.b=22
 sigma.c=222

n= a.large.number=10^5
a=rnorm(n,mu.a,sigma.a)
b=rnorm(n,mu.b,sigma.b)
c=rnorm(n,mu.c,sigma.c)
y=a+b+c
plot(density(y))
mean(y)
sd(y)

Make sure to be aware of all the assumptions we are making about y,a,b and c. If you want to do something more complex like figure out the sampling variance of the mean of y. Then do this procedure many times collecting the mean and plot it.

mysimfun=function(n,mu,sigma,stat.you.want='mean') 
   #  mu is length 3 and sigma is too.

{
n= a.large.number=10^5
    a=rnorm(n,mu[1],sigma[1])
    b=rnorm(n,mu[2],sigma[2])
    c=rnorm(n,mu[3],sigma[3])
    y=a+b+c
    plot(density(y))


return(ifelse(stat.you.want=='mean',mean(y),sd(y))
}


mu=c(mu.a,my.b,mu.c)
sigma=c(sigma.a,sigma.b,sigma.c)
mi=rep(NA,100)

Then run it in a loop of some sort.

for(i in 1:100) {mi[i]=mysimfun(10,mu,sigma,stat.you.want='mean') }

par(mfrow=c(2,1)
hist(mi)
plot(density(mi))

mean(mi)
sd(mi)

Upvotes: 3

Related Questions