I have a random variable $X$ and its pdf $f(x)=4x^3$ - which happens to be distributed as $Beta(4,1)$.

I would like to approximate $E(X)$, $Var(X)$, $P(X>0.5)$ using Monte Carlo simulation with 1,000 draws 5 times.

I know that, $$E[\bar X]=\frac{1}{n}\sum_{1=1}^n X_i=\hat \mu =\bar X \\ Var[\bar X]=\frac{1}{n-1}\sum_{1=1}^n (X_i-\bar X)^2 =\hat \sigma^2 X$$

In R, I would like to use the replicate function to simulate 1000 draws. Then

> replicate(n=1000, expr=mean(rbeta(1,shape1=4,shape2=1)), simplify = "array")

First is this the correct code? I am not summing $X_i$'s and then dividing by $n$ hoping that mean() will do both? And what is the simplify parameter (ELI5). From ETHZ,

logical or character string; should the result be simplified to a vector, matrix or higher dimensional array if possible? For sapply it must be named and not abbreviated. The default value, TRUE, returns a vector or matrix if appropriate, whereas if simplify = "array" the result may be an array of “rank” (=length(dim(.))) one higher than the result of FUN(X[[i]]).

For $Var[\bar X]$ and $P[\bar X<0.5]$ I am quite lost. Based on intuition, I will have $\bar X$ so to approximate $Var[\bar X]$, I can simulate $X_i$ via rbeta(n=1000, shape1=4, shape2=1) and putting them in a column vector $\underline{a}$. With 1000 $\bar X$ also in a column vector $\underline{b}$, I can perform $$\underline{a}-\underline{b}=\underline{c}.$$

I can square each entry of $\underline{c}$ and then sum them up and divide by $n-1$. Yet I have a suspicion that there is a 1 line solution.

Side note: I alsp know by the Central Limit Theorem, $$\bar X \dot\sim N[E(\bar X)=\mu, Var(\bar X)=\frac{\sigma^2}{n}]$$

as $n\to\infty $ or $(n\ge 30)$.


1 Answer 1


If a understand correclty you want to draw samles with 5 observations each time and calculate these estimates for every instance. First randomly draw your 5 observations

    x = rbeta(5,shape1=4,shape2=1)

The calculate the quantities you need


If you want to do this a thousand times you can indeed use

    replicate(n=1000, expr={x=rbeta(5,shape1=4,shape2=1); c(mean=mean(x), var=var(x), prob05 = mean(x<0.05))}, simplify = "array")

The simplify="array" just means that instead of a list of length 1000 with vectors of length 3 it will return a 3-by-1000 matrix


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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