R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Calculate posterior distribution for binomial model

getbinPost

R Documentation

Calculate posterior distribution for binomial model

Description

Calculate posterior distribution for the simple beta-binomial model

Usage

  getbinPost(x, object, k, n, type = c("density", "cdf"),
             rel.tol = .Machine$double.eps^0.5)

Arguments

`x`	Vector specifying, where posterior distribution should be evaluated.
`object`	An SEL object.
`k`	Number of observed successes in the binomial model.
`n`	Number of total trials.
`type`	Character specifying, whether posterior density or cdf should be evaluated.
`rel.tol`	`rel.tol` argument of the integrate subroutine, which numerically calculates the normalization constant, when a non-conjugate prior is used (ie a B-spline basis not leading to a Bernstein mixture).

Value

Numeric containing the function values corresponding to x values

Author(s)

Bjoern Bornkamp

References

Bornkamp, B. and Ickstadt, K. (2009). A Note on B-Splines for Semiparametric Elicitation. The American Statistician, 63, 373–377

Diaconis, P., and Ylvisaker, D. (1985), Quantifying Prior Opinion, Bayesian Statistics 2, (eds.) J.M. Bernardo, M.H. DeGroot, D.V. Lindley and A.F.M. Smith, Elsevier Science Publishers B.V., Amsterdam, pp. 133–156.

Examples

## Diaconis, Ylvisaker spun coin example (see references)
## simulate elicitations
x <- seq(0, 1, length=9)[2:8]
ymu <- 0.5*pbeta(x, 10, 20)+0.5*pbeta(x, 20, 10)
sig <- 0.05
set.seed(4)
y <- ymu+rnorm(7, 0, sqrt(ymu*(1-ymu))*sig)

## perform fitting with different selections
A <- SEL(x, y, d = 1, Delta = 0.001, inknts = x)
foo1 <- function(x) dbeta(x, 0.5, 0.5)
B <- SEL(x, y, d = 19, N=0, Delta = 0.02, pistar = foo1)
C <- SEL(x, y, d = 19, N=0, Delta = 0.05, pistar = foo1)
comparePlot(A, B, C, addArgs = list(layout = c(3,1)))

## posterior
sq <- seq(0,1,length=201)
res1 <- getbinPost(sq, A, 3, 10, type = "density")
res2 <- getbinPost(sq, B, 3, 10, type = "density")
res3 <- getbinPost(sq, C, 3, 10, type = "density")

## parametric posterior (corresponding to B(0.5,0.5) prior)
plot(sq, dbeta(sq, 3.5, 7.5), type="l", xlab = "",
     ylab = "Posterior density", lty = 4,
     ylim=c(0,max(c(res1, res2, res3))))
## "semiparametric" posteriors
lines(sq, res1, lty = 1)
lines(sq, res2, lty = 2)
lines(sq, res3, lty = 3)
legend(0.65,4, legend=c("Scenario A", "Scenario B", "Scenario C",
       "B(0.5,0.5) prior"), lty = 1:4)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(SEL)
Loading required package: splines
Loading required package: quadprog
Loading required package: lattice
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/SEL/getbinPost.Rd_%03d_medium.png", width=480, height=480)
> ### Name: getbinPost
> ### Title: Calculate posterior distribution for binomial model
> ### Aliases: getbinPost
> ### Keywords: misc
> 
> ### ** Examples
> 
> ## Diaconis, Ylvisaker spun coin example (see references)
> ## simulate elicitations
> x <- seq(0, 1, length=9)[2:8]
> ymu <- 0.5*pbeta(x, 10, 20)+0.5*pbeta(x, 20, 10)
> sig <- 0.05
> set.seed(4)
> y <- ymu+rnorm(7, 0, sqrt(ymu*(1-ymu))*sig)
> 
> ## perform fitting with different selections
> A <- SEL(x, y, d = 1, Delta = 0.001, inknts = x)
> foo1 <- function(x) dbeta(x, 0.5, 0.5)
> B <- SEL(x, y, d = 19, N=0, Delta = 0.02, pistar = foo1)
> C <- SEL(x, y, d = 19, N=0, Delta = 0.05, pistar = foo1)
> comparePlot(A, B, C, addArgs = list(layout = c(3,1)))
> 
> ## posterior
> sq <- seq(0,1,length=201)
> res1 <- getbinPost(sq, A, 3, 10, type = "density")
> res2 <- getbinPost(sq, B, 3, 10, type = "density")
> res3 <- getbinPost(sq, C, 3, 10, type = "density")
> 
> ## parametric posterior (corresponding to B(0.5,0.5) prior)
> plot(sq, dbeta(sq, 3.5, 7.5), type="l", xlab = "",
+      ylab = "Posterior density", lty = 4,
+      ylim=c(0,max(c(res1, res2, res3))))
> ## "semiparametric" posteriors
> lines(sq, res1, lty = 1)
> lines(sq, res2, lty = 2)
> lines(sq, res3, lty = 3)
> legend(0.65,4, legend=c("Scenario A", "Scenario B", "Scenario C",
+        "B(0.5,0.5) prior"), lty = 1:4)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>