Re: [R] ipf function in R

From: GREGOR Brian J <Brian.J.GREGOR_at_odot.state.or.us>
Date: Thu, 06 Mar 2008 08:25:51 -0800


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@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
> 

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 06 Mar 2008 - 16:29:33 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 06 Mar 2008 - 16:30:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive