# rv: a simulation-based random variable class In rv: Simulation-Based Random Variable Objects

betas <- Post$beta M <- model.matrix(fit) y.rep <- rvnorm(mean=M %**% betas, sd=sigma) mlplot(y.rep) # Summarize graphically.  Note also that sigma is also an rv object. The matrix multiplication statement returns a random vector of length \Sexpr{nrow(M)}: M %**% betas  Thus all the uncertainty in the mean estimate$X\beta$and the residual s.d. estimate$\sigma$is propagated when the replicated vector$y^\rep$is generated. In effect, this single line of code thus will in fact draw from the distribution$p({y}^\rep|y)=\int\int \Normal({y}^\rep|\mu,\sigma)p(\mu,\sigma|y)\D\mu\D\sigma.$For convenience, there is a generic method rvpredict to generate replications and predictions:  ## Replications y.rep <- rvpredict(fit)  We can also generate predictions at some other covariate values:  ## Predictions at the mean of the covariates X.pred <- data.frame(x1=mean(X$x1), x2=mean(X$x2)) y.pred <- rvpredict(fit, newdata=X.pred)  We can also perturb (add uncertainty to) the covariate x1}, then predict again.  X.rep <- X X.rep$x1 <- rnorm(n=n, mean=X.rep$x1, sd=sd(X.rep$x1))
y.pred2 <- rvpredict(fit, newdata=X.rep)
y.pred2


### Graphical summaries}

Graphical summaries are still in development, but it is now possible to plot a scatterplot with a regular vector against a random vector, showing the 50\% and 95\% uncertainty intervals along with the median, using plot(y,x,...)}, whereyis not random butxis. or we can show two random scalars plotted as a 2-dimensional scatterplot withplot(x[1],x[2],...)}.

To illustrate, let us plot the predicted intervals of the previous example, along with the data points.

Plot the predictions against y in red color; then plot the perturbed predictions with blue color.

  ## Plot predictions
plot(x = y.rep, y = D$y, rvcol="red") points(y.pred2, D$y + 0.33, rvcol="blue")


Note that the function method needs to be called explicitly to be able to plot constants vs. rv objects. If the first argument of plot(x, ...) is an rv object, one can call plot

Or, we can show a random vectors as horizontal intervals using mlplot:

  mlplot(y.rep, rvcol="red")


A histogram of the simulations of a random scalar x[1]}, can be plotted withrvhist}:

rvhist(rvnorm(mean=0, sd=1)^3, xlim=c(-3, 3),
col="red", main="Cubed standard normal")


### Example: Simulating P\'olya's Urn

This code simulates 200 iterations of the well-known P\'olya's urn problem. The parameter x/(n+1) for the Bernoulli-variate-generating function rvbern(...) is random: we can generate random variables using random parameters without much trickery; our code looks therefore more natural.

The model: \begin{eqnarray} \quad X_0 &= 1 \ \quad X_n-X_{n-1}|X_{n-1} &\sim \text{Bernoulli}( X_{n-1}/(n+1) ) \end{eqnarray}

The R code:

  x <- 1
for (n in 1:100) {
x <- x + rvbern(n=1, prob=x / (n + 1))
}

  rvhist(x / (n + 1)) # Histogram


We can see that the distribution is close to uniform, which is the limiting distribution in this case.

# Details

### Obtaining the simulation matrix

To extract the simulation matrix embedded in an rv object, use sims:

s <- sims(y.rep)
dim(s)


It is our convention to have the columns represent the random vector and the rows represent the draws from the joint distribution of the vector.

### Converting matrices and vectors of simulations to rv objects

A matrix or a vector of simulations is converted into an rv object by rvsims. Continuing the above example, we'll convert the matrix back to an rv object.

y <- rvsims(s)


You can verify that all(sims(y)==s) returns \Sexpr{all(sims(y)==s)}. % TRUE Also note that length(y) gives \Sexpr{length(y)}, since y is "just a vector."

### Coercing vectors and matrices

The function as.rv(x) coerces objects to rv objects. However, this does not mean that matrices of simulations are turned into rv objects---this is done with rvsims, as explained above. as.rv(rnorm(4000)) would return a random vector of length 4000, where each component has zero variance (and one single simulation). You probably mean rvsims(rnorm(4000)), but the correct way to generate this object is rvnorm(1).

### Obtaining simulations from R2WinBUGS

R2WinBUGS [@sturtz_2005] is an interface for calling WinBUGS within R, and obtaining the simulations as an R matrix (that is embedded in a "bugs" object). If bugsobj is the bugs object returned by the bugs(...) function call, then as.rv will coerce it into a list of random vectors, split by the parameter names: y <- as.rv(bugsobj)

### Obtaining simulations from Umacs

Umacs facilitates the construction of a Gibbs/Metropolis sampler in R \citep{Kerman:2006:umacs}, and returns the simulations wrapped in an "UmacsRun" object. Again, the coercion method as.rv will convert the Umacs object, say obj, into a list of named random vectors: y <- as.rv(obj).

# Some implementation details

rv is written in "S3" style object-oriented R rather than using the methods ("S4") package. The main reason was speed, the secondary consideration was the ease of writing new functions.

The main class is called rv. Most functions expecting an rv object have names starting with rv. For example, rvnorm, rvmean, etc.

The package also features rv-specific methods extending the basic numeric vector classes, e.g. c.rv, plot.rv, etc. However, the method dispatch in R will not call the rv class method if the first object in the argument list is not an rv object; for example, c(...) will not call c.rv in the following case: suppose that x is an object of class rv and k <- 10. Then c(k, x) will not call c.rv. To ensure the proper result, wrap the first element in as.rv:

c(as.rv(k), x) will produce a proper random vector.

# Disclaimer

This program is a work in progress, and it may contain bugs. Many new features will be eventually (and hopefully) added.

For information about random variables in R, please refer to @kerman_gelman_2007. For information about Umacs (Universal Markov chain sampler) please refer to @kerman_2006.

# References

## Try the rv package in your browser

Any scripts or data that you put into this service are public.

rv documentation built on Feb. 4, 2020, 9:10 a.m.