[R] ipf function in R

[GREGOR Brian J]

Chandra,

While I can't advise on how to use the ipf function in the cat package,
I can offer the following function that we use here to balance rebalance
arrays with new marginal totals.

#ipf.R
#Function to iteratively proportionally fit a multidimensional array
#IPF also known as Fratar method, Furness method, raking and two/three
dimensional balancing.
#This method of matrix balancing is multiplicative since the margin
factors (coefficients)
#are multiplied by the seed array to yield the balanced array. 
#Ben Stabler, benjamin.stabler at odot.state.or.us, 9.30.2003
#Brian Gregor, brian.j.gregor at odot.state.or.us, 2002

#inputs:
#1) margins_ - a list of margin values with each component equal to a
margin
#Example: row margins of 210 and 300, column margins of 140 and 370
#and third dimension margins of 170 and 340.  
#[[1]]     [,1] [,2] [,3]
#210, 300
#[[2]]
#140, 370
#[[3]]
#170, 340
#2) seedAry - a multi-dimensional array used as the seed for the IPF
#3) iteration counter (default to 100)
#4) closure criteria (default to 0.001)

#For more info on IPF see:
#Beckmann, R., Baggerly, K. and McKay, M. (1996). "Creating Synthetic
Baseline Populations."
#   Transportation Research 30A(6), 415-435.
#Inro. (1996). "Algorithms". EMME/2 User's Manual. Section 6.
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~~~~~

ipf <- function(Margins_, seedAry, maxiter=100, closure=0.001) {
    #Check to see if the sum of each margin is equal
    MarginSums. <- unlist(lapply(Margins_, sum))
    if(any(MarginSums. != MarginSums.[1])) warning("sum of each margin
not equal")

    #Replace margin values of zero with 0.001
    Margins_ <- lapply(Margins_, function(x) {
          if(any(x == 0)) warning("zeros in marginsMtx replaced with
0.001") 
          x[x == 0] <- 0.001
          x
          })

    #Check to see if number of dimensions in seed array equals the
number of
    #margins specified in the marginsMtx
    numMargins <- length(dim(seedAry))
    if(length(Margins_) != numMargins) {
        stop("number of margins in marginsMtx not equal to number of
margins in seedAry")
    }

    #Set initial values
    resultAry <- seedAry
    iter <- 0
    marginChecks <- rep(1, numMargins)
    margins <- seq(1, numMargins)

    #Iteratively proportion margins until closure or iteration criteria
are met
    while((any(marginChecks > closure)) & (iter < maxiter)) {
        for(margin in margins) {
            marginTotal <- apply(resultAry, margin, sum)
            marginCoeff <- Margins_[[margin]]/marginTotal
            marginCoeff[is.infinite(marginCoeff)] <- 0
            resultAry <- sweep(resultAry, margin, marginCoeff, "*")
            marginChecks[margin] <- sum(abs(1 - marginCoeff))
        }    
        iter <- iter + 1
    }
    
    #If IPF stopped due to number of iterations then output info
    if(iter == maxiter) cat("IPF stopped due to number of iterations\n")

    #Return balanced array
    resultAry
}


Brian Gregor, P.E.
Transportation Planning Analysis Unit
Oregon Department of Transportation
Brian.J.GREGOR at odot.state.or.us
(503) 986-4120

> 
> Message: 143
> Date: Wed, 05 Mar 2008 18:14:28 +1100
> From: Chandra Shah <Chandra.Shah at Education.monash.edu.au>
> Subject: [R] ipf function in R
> To: r-help at r-project.org
> Message-ID: <47CE4854.2020103 at Education.monash.edu.au>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
> Hi
> I have a 3 x 2 contingency table:
> 10 20
> 30 40
> 50 60
> I want to update the frequencies to new marginal totals:
> 100 130
> 40 80 110
> I want to use  the ipf (iterative proportional fitting) 
> function which 
> is apparently in the cat package.
> Can somebody please advice me how to input this data and 
> invoke ipf in R 
> to obtain an updated contingency table?
> Thanks.
> By the way I am quite new to R.
> 
> -- 
>  
> 
> Dr Chandra Shah
> Senior Research Fellow
> Monash University-ACER Centre for the Economics of Education 
> and Training
> Faculty of Education, Building 6,
> Monash University
> Victoria
> Australia 3800
> Tel. +61 3 9905 2787
> Fax +61 3 9905 9184
> 
> www.education.monash.edu.au/centres/ceet
>