# IF
## rel = 1 / (1 + w/b/k)
# AND IF 
## k ==1
# THEN
## rel = 1 / (1 + w/b)
# DEFINE
## lambda = w/b
# THEN
## rel = 1/(1+lambda)
## rel*(1+lambda) = 1
## rel+rel*lambda = 1
## rel*lambda = 1-rel
## lambda = (1-rel)/rel 
### this last formula can be used to specify what values of
### lambda should generate a nice even sampling of expected
### reliability values

library(ggplot2)
set.seed(1)

a = expand.grid(
	N = c(5,10,100)
	, lambda = (1-seq(.05,.95,.05))/seq(.05,.95,.05)
)
a$mean_rel = NA
a$sd_rel = NA
a$ci_rel_lo = NA
a$ci_rel_hi = NA
it = 1e5
pb = txtProgressBar(max=length(a[,1]),style=3)
for(i in 1:length(a[,1])){
	rels = rep(NA,it)
	for(j in 1:it){
		true_scores = rnorm(a$N[i])
		obs_1 = rnorm(a$N[i],true_scores,sqrt(a$lambda[i]))
		obs_2 = rnorm(a$N[i],true_scores,sqrt(a$lambda[i]))
		rels[j] = cor(obs_1,obs_2)
	}
	a$mean_rel[i] = mean(rels)
	a$sd_rel[i] = sd(rels)
	a$mean_rel2[i] = tanh(mean(atanh(rels)))
	a$sd_rel2[i] = tanh(sd(atanh(rels)))
	a$ci_rel_lo[i] = quantile(rels,.025)
	a$ci_rel_hi[i] = quantile(rels,.975)
	setTxtProgressBar(pb,i)
}
close(pb)

b = ddply(
	.data = a
	, .variables = .(N)
	, .fun = function(x){
		fit = lm(mean_rel~I(1/(1+lambda)),data = x)
		to_return = c(fit$coefficients,summary(fit)$r.squared)
		names(to_return) = c('int','slope','r_squared')
		return(to_return)
	}
)
print(b)

ggplot(
	data = a
	, mapping = aes(
		x = I(1/(1+lambda))
		, colour = factor(N)
		, y = mean_rel
	)
)+
geom_point()+
geom_line()