Martin,
I don't think that a doubly stochastic matrix can be obtained from an
arbitrary positive rectangular matrix. There is a theorem by Sinkhorn (Am
Math Month 1967) on the diagonal equivalence of matrices with prescribed row
and column sums. It shows that given a positive matrix A(m x n), there is a
unique matrix DAE (where D and E are m x m and n x n diagonal matrices) with
rows, k*r_i (i = 1, ..., m), and column sums, c_j (j=1,...,n) such that k =
\sum_j c_j / \sum_i r_i. Therefore, the alternative row and column
normalization algorithm (same as the iterative proportional fitting
algorithm for contingency tables) will alternate between row and column sums
being unity, while the other sum alternates between k and 1/k.
Here is a slight modification of your algorithm for the rectangular case:
bistochMat.rect <- function(m,n, tol = 1e-7, maxit = 1000) {
## Purpose: Random bistochastic *square* matrix (M_{ij}):
## M_{ij} >= 0; sum_i M_{ij} = sum_j M_{ij} = 1 (for all i,
j)
##
----------------------------------------------------------------------
## Arguments: n: (n * n) matrix dimension;
##
----------------------------------------------------------------------
## Author: Martin Maechler, Date: 16 Oct 2006, 14:47
stopifnot(maxit >= 1, tol >= 0)
M <- matrix(runif(m*n), m,n)
for(i in 1:maxit) {
rM <- rowSum(M)
M <- M / rep((rM),n)
cM <- colSum(M)
M <- M / rep((cM),each=m)
if(all(abs(c(rM,cM) - 1) < tol))
break
}
## cat("needed", i, "iterations\n")
## or rather
attr(M, "iter") <- i
M
}
Using this algorithm we get for an 8 x 4 matrix, for example, we get:
> M <- bistochMat.rect(8,4)
> apply(M,1,sum)
[1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
> apply(M,2,sum)
[1] 1 1 1 1
>
Clearly the algorithm didn't converge according to your convergence
criterion, but the row sums oscillate between 1 and 0.5, and the column sums
oscillate between 2 and 1, respectively.
It is interesting to note that the above algorithm converges if we use the
infinity norm, instead of the 1-norm, to scale the rows and columns, i.e. we
divide rows and columns by their maxima.
Best,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: [EMAIL PROTECTED]
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Martin Maechler
Sent: Monday, October 16, 2006 10:30 AM
To: Florent Bresson
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Re : Generate a random bistochastic matrix
A simpler approach --- as in similar problems ---
is to use an iterative solution which works quite fast for any 'n'.
Interestingly, the number of necessary iterations seems to
*decrease* for increasing 'n' :
bistochMat <- function(n, tol = 1e-7, maxit = 1000)
{
## Purpose: Random bistochastic *square* matrix (M_{ij}):
## M_{ij} >= 0; sum_i M_{ij} = sum_j M_{ij} = 1 (for all i,
j)
##
----------------------------------------------------------------------
## Arguments: n: (n * n) matrix dimension;
##
----------------------------------------------------------------------
## Author: Martin Maechler, Date: 16 Oct 2006, 14:47
stopifnot(maxit >= 1, tol >= 0)
M <- matrix(runif(n*n), n,n)
for(i in 1:maxit) {
M <- M / outer((rM <- rowSums(M)), (cM <- colSums(M))) * sum(rM) / n
if(all(abs(c(rM,cM) - 1) < tol))
break
}
## cat("needed", i, "iterations\n")
## or rather
attr(M, "iter") <- i
M
}
attr(M <- bistochMat(70), "iter")
## typically:
## [1] 8
attr(M <- bistochMat(10), "iter")
## -> 13, 16, 15, 17, ...
set.seed(101) ; attr(M <- bistochMat(10), "iter") ## 15
set.seed(101) ; attr(M <- bistochMat(10, tol = 1e-15), "iter")## 32
------------------------
Initially, I tried to solve the general [ n x m ] case.
It seems that needs a slightly smarter "bias correction"
instead of just 'sum(M) / n'.
If someone figures that out, please post your solution!
Regards,
Martin Maechler, ETH Zurich
______________________________________________
R-help@stat.math.ethz.ch 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.
______________________________________________
R-help@stat.math.ethz.ch 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.
