[R] Rmpfr correlation

[tri@g m@iii@g oii gvd@et@dk]

Thanks a lot - solved the issue!

BW
Troels

-----Oprindelig meddelelse-----
Fra: Rui Barradas <ruipbarradas using sapo.pt> 
Sendt: 12. juli 2020 12:59
Til: tring using gvdnet.dk; r-help mailing list <r-help using r-project.org>
Emne: Re: [R] Rmpfr correlation

Hello,

Why not write a function COR? Not one as general purpose as stats::cor but a simple one, to compute the sample Pearson correlation only.


library(Rmpfr)

COR <- function(x, y){
   precBits <- getPrec(x)[1]
   n <- mpfr(length(x), precBits = precBits)
   x.bar <- mean(x)
   y.bar <- mean(y)
   numer <- sum(x*y) - n*x.bar*y.bar
   denom <- sqrt(sum(x*x) - n*x.bar*x.bar) * sqrt(sum(y*y) - n*y.bar*y.bar)
   numer/denom
}

set.seed(2020)
KA <- mpfr(10^-4.6, 128)
x <- rnorm(100)*KA
y <- rnorm(100)*x

cor(as.numeric(x), as.numeric(y)) # -0.1874986
#[1] -0.1874986

COR(x, y)
#1 'mpfr' number of precision  128   bits
#[1] -0.1874985950531874160800643775644747505073


Hope this helps,

Rui Barradas

Às 10:42 de 12/07/20, tring using gvdnet.dk escreveu:
> Dear friends - I'm calculating buffer capacities by different methods and
> need very high precision and package Rmpfr is working beautifully. However,
> I have not been able to find out how to keep precision when finding
> correlations.
> 
> library(Rmpfr)
> 
> KA <- mpfr(10^-4.6, 128)
> 
> x <- rnorm(100)*KA
> 
> y <- rnorm(100)*x
> 
> cor(x,y) # "x" must be numeric
> 
> cor(as.numeric(x),as.numeric(y))# 0.2918954
> 
>   
> 
> In my concrete application I get cor = 1 for
> cor(as.numeric(dff$BB),as.numeric(BBVS)) even though I have
> 
>   
> 
> str(summary((dff$BB)-(BBVS)))
> Class 'summaryMpfr' [package "Rmpfr"] of length 6 and precision 128
>   4.61351010833e-8 7.33418976521e-7 1.31009046563e-5 3.76407022709e-5
> 5.72386764888e-5 ...
> 
>   
> 
> I am on windows 10
> 
> R version 3.6.1
>   
> Best wishes
> Troels Ring,
> Aalborg, Denmark
> 
>   
> 
>   
> 
> 
> This email has been scanned by BullGuard antivirus protection.
> For more info visit www.bullguard.com
> <http://www.bullguard.com/tracking.aspx?affiliate=bullguard&buyaffiliate=smt
> p&url=/>
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
> 


This email has been scanned by BullGuard antivirus protection.
For more info visit www.bullguard.com