solution6 - Homework 6 solution 1 Problem 1 Since p-value=...

Info iconThis preview shows pages 1–3. Sign up to view the full content.

View Full Document Right Arrow Icon

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Homework 6 solution 1 Problem 1 Since p-value= ∑ T ∈ Ω 1 p ( T ) ≤ p ( T ) p ( T ) = E (1 p ( T ) ≤ p ( T ) ), we can sample p ( T ) by Metropolis- Hastings algorithm, and use the sample mean of those indicator functions to approx- imate the theoretical mean. In order to estimate the se, I will estimate the mean 10 times and use the sd of those means as the estimate of the se. R code: data<-read.table("table.txt") MH1<-function(n) { x0<-data p0=1 samples<-NULL num<-0 for(i in 1:n) { k<-sample(1:12,2) l<-sample(1:12,2) x1<-x0 u1<-runif(1) if(u1<0.5) { x1[k[1],l[1]]=x1[k[1],l[1]]+1 1 x1[k[2],l[1]]=x1[k[2],l[1]]-1 x1[k[1],l[2]]=x1[k[1],l[2]]-1 x1[k[2],l[2]]=x1[k[2],l[2]]+1 } else { x1[k[1],l[1]]=x1[k[1],l[1]]-1 x1[k[2],l[1]]=x1[k[2],l[1]]+1 x1[k[1],l[2]]=x1[k[1],l[2]]+1 x1[k[2],l[2]]=x1[k[2],l[2]]-1 } if((x1[k[1],l[1]]>=0)&(x1[k[2],l[1]]>=0)&(x1[k[1],l[2]]>=0)&(x1[k[2],l[2]]>=0)) { if(u1<0.5) { r=x0[k[2],l[1]]*x0[k[1],l[2]]/(x1[k[1],l[1]]*x1[k[2],l[2]]) } else { r=x0[k[1],l[1]]*x0[k[2],l[2]]/(x1[k[2],l[1]]*x1[k[1],l[2]]) } u<-runif(1) if(u<r) { 2 x0<-x1...
View Full Document

This note was uploaded on 10/24/2010 for the course STAT 428 taught by Professor Chen during the Spring '08 term at University of Illinois, Urbana Champaign.

Page1 / 10

solution6 - Homework 6 solution 1 Problem 1 Since p-value=...

This preview shows document pages 1 - 3. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online