simTrios <- function(nfamily,pa,pmissing,RR,pedFormat=FALSE){

#print(paste("nfamily",nfamily,"pa",pa,"pmissing",pmissing,
#	"RR",RR[1],RR[2],1,NULL,sep=" "))

######### define 15 probabilities in all the 15 cells, based on probability of pa ###########

pN <- array(0,c(3,3,3)) # where subscript 3 means 0 in genotype

pN[3,3,3] = pa^0 * (1-pa)^4 * 1
pN[3,1,3] = pa^1 * (1-pa)^3 * 1
pN[3,1,1] = pa^1 * (1-pa)^3 * RR[1]
pN[1,3,3] = pa^1 * (1-pa)^3 * 1
pN[1,3,1] = pa^1 * (1-pa)^3 * RR[1]
pN[1,1,3] = pa^2 * (1-pa)^2 * 1
pN[1,1,1] = pa^2 * (1-pa)^2 * RR[1]*2
pN[1,1,2] = pa^2 * (1-pa)^2 * RR[2]
pN[3,2,1] = pa^2 * (1-pa)^2 * RR[1]
pN[2,3,1] = pa^2 * (1-pa)^2 * RR[1]
pN[1,2,1] = pa^3 * (1-pa)^1 * RR[1]
pN[1,2,2] = pa^3 * (1-pa)^1 * RR[2]
pN[2,1,1] = pa^3 * (1-pa)^1 * RR[1]
pN[2,1,2] = pa^3 * (1-pa)^1 * RR[2]
pN[2,2,2] = pa^4 * (1-pa)^0 * RR[2]

# checking
#pN[,,1]+pN[,,2]+pN[,,3]
sum(pN)

# standardization
pN <- pN / sum(pN)

# simulate multinomial distribution
N <- array(rmultinom(1,nfamily-round(nfamily*pmissing),pN),
	c(3,3,3))
M1 <- array(rmultinom(1,round(nfamily*pmissing),pN),c(3,3,3))

# define missingness for father
M=array(0,c(3,3))
for(i in 1:3){
  for (k in 1:3) {
    M[i,k]=sum(M1[i,1:3,k])
}}

trioType <- c("000","010","011","100","101","110","111","112",
	"021","201","121","122","211","212","222",
	"2?2","2?1","1?2","1?1","0?1","1?0","0?0")
MFCdata <- c(N[3,3,3], N[3,1,3], N[3,1,1], N[1,3,3], N[1,3,1],
	N[1,1,3], N[1,1,1], N[1,1,2], N[3,2,1], N[2,3,1], N[1,2,1],
	N[1,2,2], N[2,1,1], N[2,1,2], N[2,2,2],
	M[2,2], M[2,1], M[1,2], M[1,1], M[3,1], M[1,3], M[3,3])
pedData <- data.frame(trioType,MFCdata)
	
if(pedFormat==TRUE){
	pedData <- data.frame(matrix(0,nrow=3*nfamily,ncol=8))
	names(pedData) <- c("pedid","id","patid","matid","gender",
	"disease","marker1.1","marker1.2")
# Create pedigree information
	pedData$pedid <- rep(1:nfamily,each=3)
	pedData$id[3*(0:(nfamily-1))+1] <- 2*(1:nfamily)-1
	pedData$id[3*(0:(nfamily-1))+2] <- 2*(1:nfamily)
	pedData$id[3*(1:nfamily)] <- (1:nfamily)+2*nfamily
	pedData$patid[3*(1:nfamily)] <- 2*(1:nfamily)
	pedData$matid[3*(1:nfamily)] <- 2*(1:nfamily)-1
	pedData$gender <- rep(c(2,1,1,2,1,2),nfamily/2)
	pedData$disease <- rep(c(1,1,2),nfamily)
# Create markers based on MFCdata	
	pedData$marker1.1[3*(1:nfamily)-2] <- c(rep(1,MFCdata[1]+MFCdata[2]+MFCdata[3]),rep(2,MFCdata[4]+MFCdata[5]+MFCdata[6]+MFCdata[7]+MFCdata[8]),rep(1,MFCdata[9]),rep(2,MFCdata[10]+MFCdata[11]+MFCdata[12]+MFCdata[13]+MFCdata[14]+MFCdata[15]+MFCdata[16]+MFCdata[17]+MFCdata[18]+MFCdata[19]),rep(1,MFCdata[20]),rep(2,MFCdata[21]),rep(1,MFCdata[22]))
	pedData$marker1.2[3*(1:nfamily)-2] <- c(rep(1,MFCdata[1]+MFCdata[2]+MFCdata[3]+MFCdata[4]+MFCdata[5]+MFCdata[6]+MFCdata[7]+MFCdata[8]+MFCdata[9]),rep(2,MFCdata[10]),rep(1,MFCdata[11]+MFCdata[12]),rep(2,MFCdata[13]+MFCdata[14]+MFCdata[15]+MFCdata[16]+MFCdata[17]),rep(1,MFCdata[18]+MFCdata[19]+MFCdata[20]+MFCdata[21]+MFCdata[22]))
	pedData$marker1.1[3*(1:nfamily)-1] <- c(rep(1,MFCdata[1]),rep(2,MFCdata[2]+MFCdata[3]),rep(1,MFCdata[4]+MFCdata[5]),rep(2,MFCdata[6]+MFCdata[7]+MFCdata[8]+MFCdata[9]),rep(1,MFCdata[10]),rep(2,MFCdata[11]+MFCdata[12]+MFCdata[13]+MFCdata[14]+MFCdata[15]),rep(NA,MFCdata[16]+MFCdata[17]+MFCdata[18]+MFCdata[19]+MFCdata[20]+MFCdata[21]+MFCdata[22]))
	pedData$marker1.2[3*(1:nfamily)-1] <- c(rep(1,MFCdata[1]+MFCdata[2]+MFCdata[3]+MFCdata[4]+MFCdata[5]+MFCdata[6]+MFCdata[7]+MFCdata[8]),rep(2,MFCdata[9]),rep(1,MFCdata[10]),rep(2,MFCdata[11]+MFCdata[12]),rep(1,MFCdata[13]+MFCdata[14]),rep(2,MFCdata[15]),rep(NA,MFCdata[16]+MFCdata[17]+MFCdata[18]+MFCdata[19]+MFCdata[20]+MFCdata[21]+MFCdata[22]))
	pedData$marker1.1[3*(1:nfamily)] <- c(rep(1,MFCdata[1]+MFCdata[2]),rep(2,MFCdata[3]),rep(1,MFCdata[4]),rep(2,MFCdata[5]),rep(1,MFCdata[6]),rep(2,MFCdata[7]+MFCdata[8]+MFCdata[9]+MFCdata[10]+MFCdata[11]+MFCdata[12]+MFCdata[13]+MFCdata[14]+MFCdata[15]+MFCdata[16]+MFCdata[17]+MFCdata[18]+MFCdata[19]+MFCdata[20]),rep(1,MFCdata[21]+MFCdata[22]))
	pedData$marker1.2[3*(1:nfamily)] <- c(rep(1,MFCdata[1]+MFCdata[2]+MFCdata[3]+MFCdata[4]+MFCdata[5]+MFCdata[6]+MFCdata[7]),rep(2,MFCdata[8]),rep(1,MFCdata[9]+MFCdata[10]+MFCdata[11]),rep(2,MFCdata[12]),rep(1,MFCdata[13]),rep(2,MFCdata[14]+MFCdata[15]+MFCdata[16]),rep(1,MFCdata[17]),rep(2,MFCdata[18]),rep(1,MFCdata[19]+MFCdata[20]+MFCdata[21]+MFCdata[22]))
	}
list(N=nfamily,MAF=pa,PercentMissing=pmissing,
	RelativeRisk=RR,peddata=pedData)
}