inst/doc/rv-doc.R

In rv: Simulation-Based Random Variable Objects

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.pos = 'center', 
  fig.align = 'center',
  fig.height = 6,
  fig.width = 5
)

## ----echo=FALSE---------------------------------------------------------------
  set.seed(20100129)
options(digits=2)

## ----echo=TRUE----------------------------------------------------------------
library(rv)

## -----------------------------------------------------------------------------
setnsims(4000)

## -----------------------------------------------------------------------------
x <- rvnorm(mean=1:5, sd=1)

## -----------------------------------------------------------------------------
x

## -----------------------------------------------------------------------------
y <- rvpois(lambda=10)

## -----------------------------------------------------------------------------
rvmean(x)
rvsd(x)
rvquantile(x, c(0.025,0.25,0.5,0.75,0.975))
rvmedian(x)
rvmin(y)
rvmax(y)

## -----------------------------------------------------------------------------
x[3:4] <- rvbinom(size = 1, prob = c(0.1, 0.9))
x[3:4]

## -----------------------------------------------------------------------------
y <- as.rv(1:5)
y[3:4] <- x[3:4]
y

## -----------------------------------------------------------------------------
y <- (1:5)
impute(y, 3:4) <- x[3:4]
y

## -----------------------------------------------------------------------------
1/(1 + exp(-x[1]))

## -----------------------------------------------------------------------------
2*log(abs(x[2]))

## -----------------------------------------------------------------------------
x <- rvpois(lambda=1:5)
x
sort(x)
min(x)
max(x)

## -----------------------------------------------------------------------------
p <- runif(4) # Some prior probabilities.
y <- rvbinom(size=1, prob=p) # y is now a rv of length 4.
dim(y) <- c(2,2) # Make y into a 2x2 matrix.
y
y %**% y

## -----------------------------------------------------------------------------
E(y)

## -----------------------------------------------------------------------------
z <- rvnorm(1)
z > 1

## -----------------------------------------------------------------------------
Pr(z > 1)

## -----------------------------------------------------------------------------
z <- rvnorm(2)
Pr(z[1] > z[2]^2)

## -----------------------------------------------------------------------------
Pr(x[1] > x[2] & x[1] > x[4])
Pr(x[1] > x[2] | x[1] > x[4])

## -----------------------------------------------------------------------------
z <- rvnorm(n=2, mean=0, sd=1)
y <- exp(z)
y[2] <- y[2] * y[1]
x <- (y[1]-1) * (y[1]>1) * (y[2]>1)
E(x)
Pr(x>1)

## -----------------------------------------------------------------------------
  n <- 10
  ## Some covariates
  X <- data.frame(x1=rnorm(n, mean=0), x2=rpois(n, 10) - 10)
  y.mean <- (1.0 + 2.0 * X$x1 + 3.0 * X$x2)
  y <- rnorm(n, y.mean, sd=1.5) ## n random numbers
  D <- cbind(data.frame(y=y), X)
  ## Regression model fit
  fit <- lm(y ~ x1 + x2, data=D)

## -----------------------------------------------------------------------------
  Post <- posterior(fit)
  Post

## -----------------------------------------------------------------------------
sigma <- Post$sigma
betas <- Post$beta
M <- model.matrix(fit)
y.rep <- rvnorm(mean=M %**% betas, sd=sigma)
mlplot(y.rep) # Summarize graphically.

## -----------------------------------------------------------------------------
M %**% betas

## -----------------------------------------------------------------------------
  ## Replications
  y.rep <- rvpredict(fit)

## -----------------------------------------------------------------------------
  ## 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)

## -----------------------------------------------------------------------------
  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

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

## -----------------------------------------------------------------------------
  mlplot(y.rep, rvcol="red")
  mlplot(D$y, add=TRUE, col="blue", pch="x")

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

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

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

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

## -----------------------------------------------------------------------------
y <- rvsims(s)