A function to calculate David McFarland's iterative marriage model (1972)

McFarlandIterate.R

# McFarland Iterative adjustment:

McFarlandMarPredict <- function(M,uf,um,nf,nm,tol=1e-6) {
	# stick marriage mat together with those remaining unmarried
	Mmc <- rbind(cbind(M,um),c(uf,sum(uf,um)))
	# rescale rows then columns until margins add up to new margins
	for(i in 1:25){
		Mmc[-nrow(Mmc),] <- Mmc[-nrow(Mmc),]*(nm/(rowSums(Mmc)[-nrow(Mmc)]))
		Mmc[,-ncol(Mmc)] <- t(t(Mmc)[-ncol(Mmc),]*(nf/(colSums(Mmc)[-ncol(Mmc)])))
		if (sum(abs(rowSums(Mmc[-nrow(Mmc),])-nm)+abs(colSums(Mmc[,-ncol(Mmc)])-nf)) < tol) {break}
	}
	# as with M, male age rows, female age cols
	Mmc <- Mmc[-nrow(Mmc),-ncol(Mmc)]
	# male age in rows, female age in cols:
	maleRates <- Mmc/nm
	# female age in rows, male age in cols: (switched, so that all output same)
	femaleRates <- t(Mmc)/nf
	return(list(Marriages=Mmc,maleRates=maleRates,femaleRates=femaleRates))
}
# (same example data as for Henry (SE 1957 *5), from McFarland, 1972
# 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)
# need knowledge of count of intially unmarried people whose marriage are to be predicted- 
# for example just jitter data a bit
set.seed(1)
initUNMARf2 <- colSums(MAR)+unMARf+runif(4,min=-40000,max=80000)
initUNMARm2 <- rowSums(MAR)+unMARm+runif(4,min=-20000,max=90000)

McFarlandMarPredict(MAR,unMARf,unMARm,initUNMARf2,initUNMARm2)