A Function to Calculate Louis Henry's Panmictic Components for a Two-Sex Marriage Model

HenryPanmictic.R

# matrix of marriage counts with female ages in columns and male ages in rows
MAR <- matrix(c(4145,24435,8140,1865,1655,54515,45010,15030,80,6735,20870,19530,5,920,5435,42470),ncol=4)
rownames(MAR) <- colnames(MAR) 	<- c("15-19","20-24","25-29","30-60")
# unmarried males and females at start of the year
unMARf 	<- c(254876,147705,61804,415497)
unMARm 	<- c(265755,199437,114251,429655)
PanmicticRates <- function(MAR,unMARf,unMARm){
	# allocate P, array of stacked component counts
	P <- array(0,dim=c(dim(MAR),min(dim(MAR))))
	Pi <- Ri <- MAR
	
	# scale residuals by proportions of submatrices going iteratively down the diagonal
	for (i in 1:(min(dim(MAR))-1)){
		P[,,i] 	<-	Pi
		P[(i+1):nrow(P),(i+1):ncol(P),i] <- Pi[(i+1):nrow(P),i]%o%Pi[i,(i+1):ncol(P)]/Pi[i,i]
		Pi <- Ri <- Ri - P[,,i]
	}
	
	# component counts, stacked
	P[,,min(dim(MAR))] <- Pi
	
	# component rates
	Mcomp <- apply(P,3,rowSums)/unMARm
	Fcomp <- apply(P,3,colSums)/unMARf
	rownames(Mcomp) <- rownames(Fcomp) <- rownames(MAR)
	colnames(Mcomp) <- colnames(Fcomp) <- paste(1:ncol(Mcomp),"comp",sep="")
	
	return(list(Mcomp=Mcomp,Fcomp=Fcomp))
}
# remember females columns, males rows
PanmicticRates(MAR,unMARf,unMARm)