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, tr...@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@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.
______________________________________________
R-help@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.
