如何将Rmpfr值放入R中的函数中? [英] How to put Rmpfr values into a function in R?
问题描述
我正在计算Vandermonde矩阵的逆.我已经编写了如下代码来通过其公式显式计算逆:
I am calculating the inverse of a Vandermonde Matrix. I have written the codes to calculate the inverse explicitly by its formula as below:
library(gtools)
#input is the generation vector of terms of Vandermonde matrix.
FMinv <- function(base){
n=length(base)
inv=matrix(nrow=n,ncol=n)
for (i in 1:n){
for (j in 1:n){
if(j<n){
a=as.matrix(combinations(n,n-j,repeats.allowed = F))
arow.tmp=nrow(a) #this is in fact a[,1]
b=which(a==i)%%length(a[,1])
nrowdel=length(b)
b=replace(b,b==0,length(a[,1]))
a=a[-b,]
if(arow.tmp-nrowdel>1){
a=as.matrix(a)
nrowa=nrow(a)
prod=vector()
for(k in 1:nrowa){
prod[k]=prod(base[a[k,]])
}
num=sum(prod)
}
if(arow.tmp-nrowdel==1){
num=prod(base[a])
}
den=base[i]*prod(base[-i]-base[i])
inv[i,j]=(-1)^(j-1)*num/den
}
if(j==n){
inv[i,j]=1/(base[i]*prod(base[i]-base[-i]))
}
}
}
return(inv)
}
我定义一个基数如下:
> library(Rmpfr)
> a=mpfr(c(10:1),1000)/Rmpfr::mpfr(sum(1:10),1000)
> a
10 'mpfr' numbers of precision 1000 bits
[1] 0.18181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181819
[2] 0.16363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363637
[3] 0.14545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545456
[4] 0.12727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727274
[5] 0.10909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909091
[6] 0.090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909090909094
[7] 0.072727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727272727278
[8] 0.054545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545455
[9] 0.036363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363636363639
[10] 0.018181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181818181819
但是,当我尝试将"a"放入函数中时,我得到了:
However, when I was attempting to put the "a" into the function, I got:
> FMinv(a)
Error in sum(prod) : invalid 'type' (list) of argument
通过检查其类型,
> typeof(a)
[1] "list"
我唯一知道将其转换为值的是Rmpfr中的asNumeric().但是,
The only thing that I know to transform it to values is asNumeric() in Rmpfr. However,
> asNumeric(a)
[1] 0.18181818 0.16363636 0.14545455 0.12727273 0.10909091 0.09090909 0.07272727 0.05454545 0.03636364 0.01818182
它丢失了剩余的数字.
是否有将"a"放入函数中而不会丢失小数?
Is there anyway to put the "a" into my function without losing decimals?
谢谢!
推荐答案
诀窍是使用S3方法.
定义一个通用的默认方法,该默认方法将使用您的正常"数字(即类"numeric"
的对象)和问题所要查询的函数来调用.
The trick is to use S3 methods.
Define a generic, a default method to be called with your "normal" numbers, meaning, objects of class "numeric"
and the function the question is asking for.
那是问题之一.花了一段时间,但我相信下面的代码是正确的.
That is the problem one. It took a while but I believe the code below is right.
library(OBsMD)
FMinv <- function(...) UseMethod("FMinv")
FMinv.default <- function(base) {
# Your function
# unchanged
}
FMinv.mpfr <- function(base, precBits = getPrec(base)) {
n <- length(base)
inv <- mpfr(rep(0, n*n), precBits = precBits)
inv <- matrix(inv, nrow = n, ncol = n)
for (i in 1:n) {
for (j in 1:n) {
if (j < n) {
a <- combinations(n, n - j, repeats.allowed = F)
a <- as.matrix(a)
arow.tmp <- nrow(a) # this is in fact a[, 1]
b <- which(a == i) %% length(a[, 1])
nrowdel <- length(b)
b <- replace(b, b == 0, length(a[, 1]))
a <- a[-b, ]
num <- mpfr(0, precBits[1])
if (arow.tmp - nrowdel > 1) {
a <- as.matrix(a)
nrowa <- nrow(a)
for (k in 1:nrowa) {
num <- num + prod(base[a[k, ]])
}
}
if (arow.tmp - nrowdel == 1) {
num <- num + prod(base[a])
}
den <- base[i] * prod(base[-i] - base[i])
inv[i, j] <- (-1)^(j - 1) * num/den
}
if (j == n) {
inv[i, j] <- 1/(base[i] * prod(base[i] - base[-i]))
}
}
}
return(inv)
}
现在测试这两种方法并比较一些结果的值.
Now test both methods and compare some of the results' values.
library(Rmpfr)
a <- mpfr(c(10:1),1000)/Rmpfr::mpfr(sum(1:10),1000)
inv1 <- FMinv(asNumeric(a))
inv2 <- FMinv(a)
inv1[10, 10]
#[1] -6.98014e+11
inv2[10, 10]
#1 'mpfr' number of precision 1000 bits
#[1] -698013564040.84166942239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906525573192239858906474
这篇关于如何将Rmpfr值放入R中的函数中?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!