Hi, Ravi:
Question: Are duplicate real numbers complex conjugates? For
some purposes, perhaps. However, for identifying (damped) periodicity
on a difference or differential equation, we must exclude repeated real
roots.
In any event, your suggestion inspired me to create a function
"findConjugates" (see below), which I've added to the "FinTS" package
currently under development on R-Forge. The source code is currently
available via "svn checkout
svn://svn.r-forge.r-project.org/svnroot/fints". In another day or so,
it should be available (with documentation) via
'install.packages("FinTS",repos="http://r-forge.r-project.org";)'. In
another couple of months, it should appear on CRAN.
Thanks again for your suggestion.
Best Wishes,
Spencer
######################################
findConjugates <- function(x,
complex.eps=1000*.Machine[["double.neg.eps"]]){
##
## 1. compute normalization
##
if(length(x)<1)return(complex(0))
ax <- abs(x)
m2 <- outer(ax, ax, pmax)
##
## 2. Compute complex differences
##
c2 <- (abs(outer(x, Conj(x), "-") / m2) < complex.eps)
c2[m2==0] <- FALSE
c2 <- (c2 & lower.tri(c2))
##
## 3. Any differences exceed complex.eps?
##
if(any(c2)){
# check standard differences
d2 <- (abs(outer(x, x, "-") / m2) > complex.eps)
d2[m2==0] <- FALSE
#
cd2 <- (c2 & d2)
if(any(cd2)){
ic <- sort(unique(row(cd2)[cd2]))
return(x[ic])
}
}
complex(0)
}
######################################
Ravi Varadhan wrote:
> Hi Spencer,
>
> Here is a simple approach to detect conjugate pairs:
>
> is.conj <- function(z1, z2, tol=1.e-10) {
> # determine if two complex numbers are conjugates
> cond1 <- abs(Re(z1) - Re(z2)) < tol
> cond2 <- abs(Im(z1) + Im(z2)) < tol
> cond1 & cond2
> }
>
> set.seed(123)
> z <- polyroot(sample(1:5, size=8, rep=T))
> zmat <- which(outer(z, z, FUN="is.conj"), arr.ind=T)
> zmat[zmat[,1] < zmat[,2], ]
>
> # result
> row col
> [1,] 1 3
> [2,] 5 6
> [3,] 4 7
>
>
> We see that (1,3), (4,7), and (5,6) are the conjugate pairs.
>
> This doesn't address the issue of numerical round-off (there is no argument
> in polyroot that governs the accuracy of the roots).
>
> 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 Spencer Graves
> Sent: Friday, November 23, 2007 12:08 PM
> To: r-help@r-project.org
> Subject: [R] complex conjugates roots from polyroot?
>
> Hi, All:
>
> Is there a simple way to detect complex conjugates in the roots
> returned by 'polyroot'? The obvious comparison of each root with the
> complex conjugate of the next sometimes produces roundoff error, and I
> don't know how to bound its magnitude:
>
> (tst <- polyroot(c(1, -.6, .4)))
> tst[-1]-Conj(tst[-2])
> [1] 3.108624e-15+2.22045e-16i
> abs(tst[-1]-Conj(tst[-2]))/abs(tst[-1])
> 1.971076e-15
> .Machine$double.neg.eps
> 1.110223e-16
>
> Testing (abs(tst[-1]-Conj(tst[-2]))/abs(tst[-1]) <
> (20*.Machine$double.neg.eps)) would catch this example, but it might not
> catch others.
>
> The 'polyroot' help page says, "See Also ... the 'zero' example in
> the demos directory." Unfortunately, I've so far been unable to find
> that.
>
> Any suggestions?
> Thanks in advance.
> Spencer Graves
>
> ______________________________________________
> R-help@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.
>
> ______________________________________________
> R-help@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.
>
______________________________________________
R-help@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.
