[R] complex conjugates roots from polyroot?
Spencer Graves
spencer.graves at pdf.com
Mon Nov 26 20:01:21 CET 2007
Hi, Ravi:
Thanks very much. Per your suggestion, I made
complex.eps=.Machine[["double.eps"]]), added your examples to the help
file. I also listed "author" as "Spencer Graves and Ravi Varadhan" in
the help file and uploaded the changes to R-forge, as I mentioned below.
Later, I will modify the code that computes the roots to use
roots(polynomial(...)), per your suggestion.
Best Wishes,
Spencer
Ravi Varadhan wrote:
> Spencer,
>
> I just observed that the polynomial root calculation in the package,
> "polynom", using the function solve() is more accurate than the polyroot()
> function in the "base" package. Here is an example:
>
> set.seed(1234)
> p <- polynomial(sample(1:10, size=45, rep=T)) # degree 44
> z <- solve(p)
> findConjugates(z, complex.eps=.Machine$double.eps) # this identifies all 21
> conjugate pairs
>
> z1 <- polyroot(p)
> findConjugates(z1, complex.eps=.Machine$double.eps) # this only identifies
> only 3 conjugate pairs
>
> As, I had mentioned earlier, I can't tell what tolerances are used by these
> algorithms. However, it appears that "polynom" is more accurate.
>
> 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: rvaradhan at jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----------------------------------------------------------------------------
> --------
>
> -----Original Message-----
> From: Spencer Graves [mailto:spencer.graves at pdf.com]
> Sent: Sunday, November 25, 2007 12:36 PM
> To: Ravi Varadhan
> Cc: r-help at r-project.org
> Subject: Re: [R] complex conjugates roots from polyroot?
>
> 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: rvaradhan at jhmi.edu
>>
>> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>>
>>
>>
>>
>>
> ----------------------------------------------------------------------------
>
>> --------
>>
>> -----Original Message-----
>> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
>>
> On
>
>> Behalf Of Spencer Graves
>> Sent: Friday, November 23, 2007 12:08 PM
>> To: r-help at 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 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.
>>
>> ______________________________________________
>> 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.
>>
>>
>
>
>
More information about the R-help
mailing list