richarddmorey

Nmin = 20
Nmax = 100
Mmax = 100
Mtotal = 2*Nmin*Mmax

Num_sim_per_N = 10

N = rep(seq(Nmin,Nmax,10), each = Num_sim_per_N)
M = round( Mtotal / ( 2 * N ) )

richarddmorey

## This function computes the (normalized) marginal likelihood.
## You will not use this function directly.
p.y.alt = function(t.stat, N, log.prior.dens,lo=-Inf,up=Inf,...){
  normalize = integrate(function(delta,...){
    exp(log.prior.dens(delta, ...))
  }, lower = lo, upper = up, ...)[[1]]
  py = integrate(function(delta,t.stat,N,...){
    exp(
      dt(t.stat, N - 1, ncp = delta*sqrt(N), log = TRUE) +

richarddmorey

x = scan()
35.84 8.71 38.14 8.55 -3.89 36 .000
38.22 7.66 39.73 7.80 -5.19 36 .000
67.16 8.05 68.65 8.78 -4.08 36 .000
22.08 3.28 22.68 3.58 -2.12 36 .041
22.30 2.46 23.03 2.65 -3.96 36 .000
22.78 2.97 22.95 2.91 -0.92 36 .362
60.73 2.75 63.76 2.80 -5.85 36 .000
18.19 1.68 19.84 1.52 -4.84 36 .000

richarddmorey

############
# Option 1: Get all solutions
# using brute force
############

## This function creates every possible distribution of responses for
## a likert scale with nlev responses. This is total brute force. There's
## probably a better way.
## Argument:
## v : initially, the total number of responses

richarddmorey

###########
# S4 class to represent real numbers as logarithms
# Includes standard arithmetic operations +,-,*,/, and ^
# Richard D. Morey (richarddmorey@gmail.com)
# November 2014
###########

# Compute log(exp(a) + exp(b))
logExpAplusExpB <- function(a, b)
{

richarddmorey

# Install the BayesFactor and devtools packages, if you don't already have them

# Load my S4 class for representing real numbers with logarithms
# and performing arithmetic on them
# See the code at https://gist.github.com/richarddmorey/3c77d0065983e31241bff3807482443e
devtools::source_gist('3c77d0065983e31241bff3807482443e')

# set random seed so results are reproducible
set.seed(2)

richarddmorey

\documentclass{article}
\usepackage{fancyvrb, listings, xcolor}
\usepackage{hyperref} 

\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{mauve}{rgb}{0.88, 0.69, 1.0}

\lstset{ %
  language=R,                     % the language of the code
  basicstyle=\footnotesize,       % the size of the fonts that are used for the code

richarddmorey

---
title: "Joint and conditional, binomial/beta example"
author: "Richard D. Morey"
date: "4 June 2016"
output: html_document
---

```{r,echo=FALSE,warning=FALSE}
library(rgl)
library(knitr)

richarddmorey

devtools::source_gist("e49c345fcd6cfe8f32535eda4bd8c29b", filename='utility_functions.R')

# Prior Setup
p0 = .5
rscale = .5
interval = c(.5,1)
# Do not change shift unless you want different prior shift
# will cause deviation from BayesFactor package
# (not that there's anything wrong with that, but changes demo)
shift = qlogis(p0)

richarddmorey

y = 13
N = 35

a = 1
b = 1

theta = seq(0, 1, len = 100)
likelihood = theta ^ y * (1-theta)^(N-y) #dbinom(y, N, theta)
prior = theta ^ (a-1) * (1-theta)^(b-1) #dbeta(theta, a, b)
posterior1 = prior * likelihood