#

####################################
## simulating the Singles method
## for estimating segregation ratio
####################################

p.hat_rep(0,10)
p_rep(0,10)

# we will repeat the whole process 10 times to get
# and average and a sample sd
for(k in 1:10){

# create a "population" of 5 child families
# whose parents can have affected childeren with
# probability 1/4
mat <- matrix(rbinom(50000,1,.25),ncol=5)

# create a proband vector = # of probands in each family where
# the ascertainment prob for individual = .05
proband <- rep(0,10000)
for(i in 1:10000){
if(sum(mat[i,])>0){
proband[i] <- rbinom(1,sum(mat[i,]),.05)
}#if
}#for

# create matrix of selected families (at least one proband)
selected.fams <- mat[proband>=1,]
proband <- proband[proband>=1]

d_sum(proband==1)

# vector of number affected in selected.fams

num.aff <- rep(0,dim(selected.fams)[1])
for(i in 1:dim(selected.fams)[1]){
num.aff[i] <- sum(selected.fams[i,])
}#for

# find out estimate for p
num <- (sum(num.aff)-d)
denom <- (dim(selected.fams)[1]*dim(selected.fams)[2]-d)

p.hat[k] <- (num/denom)

# what was actual p?  should be about .25

p[k] <- (sum(mat)/(dim(mat)[1]*dim(mat)[2]))

}#for

### OUTPUT ##################################################################

> p.hat

[1] 0.2564184 0.2508774 0.2549020 0.2450067 0.2455657 0.2479348 0.2508532

[8] 0.2550255 0.2562375 0.2509689

> p

[1] 0.25202 0.24916 0.24856 0.24950 0.24934 0.24774 0.24784 0.25038 0.25008

[10] 0.24730

> p.hat-p

[1]  0.0043983502  0.0017173796  0.0063419608 -0.0044933422 -0.0037742845

[6]  0.0001947709  0.0030132423  0.0046455352  0.0061574987  0.0036688908

> stdev(p.hat-p)

[1] 0.003818014

> mean(p.hat-p)

[1] 0.002187

##############################################################################