Homework 4
Kate Lindsey
January 31, 2015
Sara Kessler and I struggled for twelve hours to complete this assignment. Our solution to problem 1 is completely wrong, and we could not figure out why. We know it’s wrong because our confidence intervals are getting larger, not smaller, and they also don’t center around .4 or .6. Our solution to problem 2 is better. Our solution to problem 3 also at least looks like it’s giving us reasonable data. I really hope we go over this material in class!
1
table=matrix(0,20,7)
colnames(table)=c('n','LCI1','UCI1','LCI2','UCI2','overlap','bayes.factor')
parameters = seq(from=0, to=1, by=.01)
parameters## [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13
## [15] 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27
## [29] 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41
## [43] 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55
## [57] 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69
## [71] 0.70 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83
## [85] 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97
## [99] 0.98 0.99 1.00prior = rep(1/length(parameters), length(parameters))
prior## [1] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [7] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [13] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [19] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [25] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [31] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [37] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [43] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [49] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [55] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [61] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [67] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [73] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [79] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [85] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [91] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099
## [97] 0.00990099 0.00990099 0.00990099 0.00990099 0.00990099num.observed=seq(5,100,5)
num.observed## [1] 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85
## [18] 90 95 100for(i in 1:20) {
num.passive.A = sum(rbinom(num.observed[i],1,.4))
likelihood.A=sapply(parameters,function(p) dbinom(num.passive.A, size=num.observed[i],prob=p))
non.normalized.posterior.A=likelihood.A*prior
posterior.A = non.normalized.posterior.A/sum(non.normalized.posterior.A)
lower.CI.A=quantile(posterior.A,0.025)
upper.CI.A=quantile(posterior.A,0.975)
num.passive.B = sum(rbinom(num.observed[i],1,.6))
likelihood.B=sapply(parameters,function(p) dbinom(num.passive.B, size=num.observed[i],prob=p))
non.normalized.posterior.B=likelihood.B*prior
posterior.B = non.normalized.posterior.B/sum(non.normalized.posterior.B)
lower.CI.B=quantile(posterior.B,0.025)
upper.CI.B=quantile(posterior.B,0.975)
overlap=(upper.CI.A-lower.CI.B>0)
bayes.factor=((dbinom(num.passive.A,size=num.observed[i],.5)*dbinom(num.passive.B,size=num.observed[i],.5))/(dbinom(num.passive.A,size=num.observed[i],.4)*dbinom(num.passive.B,size=num.observed[i],.6)))
table[i,]=c(num.observed[i],lower.CI.A,upper.CI.A,lower.CI.B,upper.CI.B,overlap,bayes.factor)
}
table## n LCI1 UCI1 LCI2 UCI2 overlap
## [1,] 5 2.598990e-06 0.02068247 2.501125e-08 0.02448941 1
## [2,] 10 2.767661e-16 0.04210868 2.636144e-09 0.02693533 1
## [3,] 15 1.819789e-13 0.03280051 4.798272e-12 0.03226423 1
## [4,] 20 3.012951e-15 0.03680966 3.508890e-16 0.03663362 1
## [5,] 25 1.180340e-35 0.06036774 1.451203e-22 0.04217187 1
## [6,] 30 1.871094e-22 0.04455510 3.877481e-30 0.04791375 1
## [7,] 35 2.773186e-24 0.04760327 9.887275e-32 0.04983532 1
## [8,] 40 6.776738e-28 0.05063029 1.564627e-31 0.05137380 1
## [9,] 45 5.779209e-39 0.05554977 6.819228e-31 0.05341724 1
## [10,] 50 2.100333e-33 0.05602429 1.134290e-40 0.05743687 1
## [11,] 55 2.006088e-42 0.05948886 8.163293e-40 0.05841585 1
## [12,] 60 7.608536e-50 0.06300741 6.159606e-48 0.06250135 1
## [13,] 65 7.966267e-45 0.06322267 1.426938e-51 0.06481735 1
## [14,] 70 3.295321e-55 0.06701059 2.696776e-49 0.06539507 1
## [15,] 75 2.732758e-53 0.06811861 5.702371e-57 0.06852177 1
## [16,] 80 6.317057e-57 0.07021347 3.014973e-68 0.07229241 1
## [17,] 85 8.753709e-59 0.07138597 3.453928e-55 0.07147214 1
## [18,] 90 6.852555e-68 0.07437072 9.153691e-70 0.07477892 1
## [19,] 95 7.202631e-62 0.07524303 2.677259e-75 0.07692890 1
## [20,] 100 1.760679e-71 0.07741131 1.005992e-84 0.07982583 1
## bayes.factor
## [1,] 0.545081342
## [2,] 0.297113670
## [3,] 1.229816301
## [4,] 1.508287321
## [5,] 0.014257153
## [6,] 0.672199171
## [7,] 0.824407259
## [8,] 0.449369016
## [9,] 0.244942666
## [10,] 0.300405774
## [11,] 0.828962011
## [12,] 0.026445826
## [13,] 0.369443918
## [14,] 0.453098220
## [15,] 0.109766838
## [16,] 0.017727957
## [17,] 4.231446910
## [18,] 0.017776862
## [19,] 0.110373287
## [20,] 0.005281749plot(table[,7],num.observed,xlab="bayes.factor")
conf.int.width.A=table[,3]-table[,2]
conf.int.width.B=table[,5]-table[,4]
plot(conf.int.width.A,num.observed)
plot(conf.int.width.B,num.observed)
2
table.rej=matrix(0,20,7)
colnames(table.rej)=c('n','LCI1','UCI1','LCI2','UCI2','overlap','bayes.factor')
n.observations = seq(5,100,5)
for(i in 1:20) {
accepted.parameters.A = c()
accepted.parameters.B = c()
passives.in.observation.A = sum(rbinom(n=1,size=n.observations[i],prob=.4))
passives.in.observation.B = sum(rbinom(n=1,size=n.observations[i],prob=.6))
while (length(accepted.parameters.A) < 500) {
# make sure to get 500 samples from the conditional distribution
this.sim.param.A = runif(1,0,1)
passives.in.simulated.data.A = rbinom(n=1, size=n.observations[i], prob=this.sim.param.A)
if (passives.in.simulated.data.A == passives.in.observation.A) {
accepted.parameters.A = c(accepted.parameters.A, this.sim.param.A)
}
}
lower.CI.A = quantile(accepted.parameters.A, .025)
upper.CI.A = quantile(accepted.parameters.A, .975)
while (length(accepted.parameters.B) < 500) {
# make sure to get 500 samples from the conditional distribution
this.sim.param.B = runif(1,0,1)
passives.in.simulated.data.B = rbinom(n=1, size=n.observations[i], prob=this.sim.param.B)
if (passives.in.simulated.data.B == passives.in.observation.B) {
accepted.parameters.B = c(accepted.parameters.B, this.sim.param.B)
}
}
lower.CI.B = quantile(accepted.parameters.B, .025)
upper.CI.B = quantile(accepted.parameters.B, .975)
overlap=(upper.CI.A-lower.CI.B>0)
bayes.factor=((dbinom(passives.in.observation.A,size=n.observations[i],.5)*dbinom(passives.in.observation.B,size=n.observations[i],.5))/(dbinom(passives.in.observation.A,size=n.observations[i],.4)*dbinom(passives.in.observation.B,size=n.observations[i],.6)))
table.rej[i,]=c(n.observations[i],lower.CI.A,upper.CI.A,lower.CI.B,upper.CI.B,overlap,bayes.factor)
}
table.rej## n LCI1 UCI1 LCI2 UCI2 overlap bayes.factor
## [1,] 5 0.2334490 0.8660126 0.2388867 0.8941807 1 1.226433020
## [2,] 10 0.3110727 0.8315390 0.2521010 0.7576313 1 2.256206929
## [3,] 15 0.1131400 0.5183364 0.4369743 0.8375225 1 0.161951118
## [4,] 20 0.1544866 0.5253775 0.4266386 0.8106774 1 0.132414799
## [5,] 25 0.2609896 0.6306326 0.4016905 0.7685852 1 0.548092851
## [6,] 30 0.2092982 0.5448136 0.4894991 0.8074435 1 0.088520055
## [7,] 35 0.2299680 0.5386842 0.4068043 0.7090235 1 0.244268818
## [8,] 40 0.2657982 0.5651944 0.4281418 0.7166839 1 0.299579344
## [9,] 45 0.2529587 0.5148502 0.5562099 0.8277219 0 0.014335922
## [10,] 50 0.2555457 0.5093820 0.5221098 0.7701130 0 0.026373072
## [11,] 55 0.2306471 0.4759970 0.3404193 0.6020839 1 0.552641341
## [12,] 60 0.3436590 0.5957578 0.5725033 0.8021781 1 0.039668739
## [13,] 65 0.1961556 0.4080208 0.5120804 0.7430245 0 0.001898288
## [14,] 70 0.2552557 0.4797929 0.4226265 0.6458146 1 0.089500883
## [15,] 75 0.3487848 0.5589505 0.5175094 0.7240422 1 0.109766838
## [16,] 80 0.2894615 0.4863854 0.4557829 0.6622235 1 0.089747783
## [17,] 85 0.3485093 0.5476023 0.5121394 0.7114724 1 0.073379763
## [18,] 90 0.3091711 0.4991984 0.5357340 0.7322791 0 0.005267218
## [19,] 95 0.3283880 0.5178126 0.4853033 0.6921642 1 0.049054794
## [20,] 100 0.3715885 0.5328705 0.5235969 0.7034089 1 0.060162419plot(table.rej[,7],n.observations, xlab="bayes.factor")
conf.int.width.A=table.rej[,3]-table.rej[,2]
conf.int.width.B=table.rej[,5]-table.rej[,4]
plot(conf.int.width.A,n.observations)
plot(conf.int.width.B,n.observations)
3
table.nup=matrix(0,20,7)
colnames(table.nup)=c('n','LCI1','UCI1','LCI2','UCI2','overlap','bayes.factor')
n.observations.nup = seq(5,100,5)
for(i in 1:20) {
accepted.parameters.A.nup = c()
accepted.parameters.B.nup = c()
sampler = function(j) {
q = runif(1,0,1)
if (q < .7) {
parameter = rbinom(1,1,.95)
} else {
if (q < .72) {
parameter = rbinom(1,1,1)
} else {
if (q < .74) {
parameter = rbinom(1,1,0)
} else {
parameter = rbinom(1,1,.5)
}}}
simulated.heads=rbinom(n=1, size=n.observations.nup[i],prob=parameter)
return(sum(simulated.heads))
}
obs.set.A.nup = sampler()
obs.set.B.nup = sampler()
while (length(accepted.parameters.A.nup) < 500) {
# make sure to get 500 samples from the conditional distribution
this.sim.param.A.nup = runif(1,0,1)
simulated.flips.A.nup = rbinom(n=1, size=n.observations.nup[i], prob=this.sim.param.A.nup)
if (simulated.flips.A.nup == obs.set.A.nup) {
accepted.parameters.A.nup = c(accepted.parameters.A.nup, this.sim.param.A.nup)
}
}
lower.CI.A.nup = quantile(accepted.parameters.A.nup, .025)
upper.CI.A.nup = quantile(accepted.parameters.A.nup, .975)
while (length(accepted.parameters.B.nup) < 500) {
# make sure to get 500 samples from the conditional distribution
this.sim.param.B.nup = runif(1,0,1)
simulated.flips.B.nup = rbinom(n=1, size=n.observations.nup[i], prob=this.sim.param.B.nup)
if (simulated.flips.B.nup == obs.set.B.nup) {
accepted.parameters.B.nup = c(accepted.parameters.B.nup, this.sim.param.B.nup)
}
}
lower.CI.B.nup = quantile(accepted.parameters.B.nup, .025)
upper.CI.B.nup = quantile(accepted.parameters.B.nup, .975)
overlap.nup=(upper.CI.A.nup-lower.CI.B.nup>0)
bayes.factor.nup=((dbinom(obs.set.A.nup,size=n.observations.nup[i],.5)*dbinom(obs.set.B.nup,size=n.observations.nup[i],.5))/(dbinom(obs.set.A.nup,size=n.observations.nup[i],.4)*dbinom(obs.set.B.nup,size=n.observations.nup[i],.6)))
table.nup[i,]=c(n.observations.nup[i],lower.CI.A.nup,upper.CI.A.nup,lower.CI.B.nup,upper.CI.B.nup,overlap.nup,bayes.factor.nup)
}
table.nup## n LCI1 UCI1 LCI2 UCI2 overlap
## [1,] 5 0.5953192678 0.99389939 0.5257521102 0.99403327 1
## [2,] 10 0.7112735078 0.99670186 0.7240572347 0.99796034 1
## [3,] 15 0.7981326441 0.99826843 0.8004282698 0.99894903 1
## [4,] 20 0.8383382309 0.99910151 0.8508520968 0.99883730 1
## [5,] 25 0.8729308653 0.99915769 0.0014779890 0.12205553 1
## [6,] 30 0.8914601841 0.99881007 0.8802733059 0.99905409 1
## [7,] 35 0.8939022869 0.99867285 0.9026553892 0.99908435 1
## [8,] 40 0.9207051740 0.99933398 0.9224680258 0.99958598 1
## [9,] 45 0.9324904339 0.99965118 0.9254681586 0.99937941 1
## [10,] 50 0.9376058638 0.99919640 0.9251270799 0.99975933 1
## [11,] 55 0.9359354024 0.99961649 0.9344083190 0.99946740 1
## [12,] 60 0.0003968360 0.05899606 0.0004182535 0.05930944 1
## [13,] 65 0.9473123594 0.99949413 0.9505104445 0.99960366 1
## [14,] 70 0.9493512598 0.99976217 0.9502276835 0.99942578 1
## [15,] 75 0.9537807362 0.99959220 0.9546617364 0.99977802 1
## [16,] 80 0.9495842992 0.99953720 0.0003405729 0.04665579 1
## [17,] 85 0.0003823705 0.04293921 0.0002403372 0.04290507 1
## [18,] 90 0.9656385927 0.99971708 0.9604061951 0.99955476 1
## [19,] 95 0.0004957268 0.04025073 0.0003366174 0.03741153 1
## [20,] 100 0.9675069391 0.99973880 0.9626862145 0.99976764 1
## bayes.factor
## [1,] 1.226433e+00
## [2,] 1.504138e+00
## [3,] 1.844724e+00
## [4,] 2.262431e+00
## [5,] 7.006492e+04
## [6,] 3.403008e+00
## [7,] 4.173562e+00
## [8,] 5.118594e+00
## [9,] 6.277613e+00
## [10,] 7.699071e+00
## [11,] 9.442395e+00
## [12,] 1.158047e+01
## [13,] 1.420267e+01
## [14,] 1.741862e+01
## [15,] 2.136277e+01
## [16,] 3.203333e+15
## [17,] 3.213255e+01
## [18,] 3.940842e+01
## [19,] 4.833179e+01
## [20,] 5.927570e+01plot(table.nup[,7],n.observations.nup,xlab="bayes.factor")
conf.int.width.A.nup=table.nup[,3]-table.nup[,2]
conf.int.width.B.nup=table.nup[,5]-table.nup[,4]
plot(conf.int.width.A.nup,n.observations.nup)
plot(conf.int.width.B.nup,n.observations.nup)