Re: [Rd] Discrepancy: R sum() VS C or Fortran sum
My simple functions were to compare the result with the gfortran
compiler sum() function. I thought that the Fortran sum could not be
less precise than R. I was wrong. I am impressed. The R sum does in fact
match the result if we use the Kahan algorithm.
P.
I am glad to see that R sum() is more accurate than the gfortran
compiler sum.
On 16/03/18 11:37 AM, luke-tier...@uiowa.edu wrote:
Install the gmp package, run your code, and then try this:
bu <- gmp::as.bigq(u)
bs4 <- bu[1] + bu[2] + bu[3] + bu[4] + bu[5]
s4 <- as.double(bs4)
s1 - s4
## [1] 0
s2[[2]] - s4
## [1] 7.105427e-15
s3 - s4
## [1] 7.105427e-15
identical(s1, s4)
## [1] TRUE
`bs4` is the exact sum of the binary rationals in your `u` vector;
`s4` is the closest double precision to this exact sum.
Looking at the C source code for sum() will show you that it makes
some extra efforts to get a more accurate sum than your simple
version.
Best,
luke
On Fri, 16 Mar 2018, Pierre Chausse wrote:
Hi all,
I found a discrepancy between the sum() in R and either a sum done in
C or Fortran for vector of just 5 elements. The difference is very
small, but this is a very small part of a much larger numerical
problem in which first and second derivatives are computed
numerically. This is part of a numerical method course I am teaching
in which I want to compare speeds of R versus Fortran (We solve a
general equilibrium problem all numerically, if you want to know).
Because of this discrepancy, the Jacobian and Hessian in R versus in
Fortran are quite different, which results in the Newton method
producing a different solution (for a given stopping rule). Since the
solution computed in Fortran is almost identical to the analytical
solution, I suspect that the sum in Fortran may be more accurate
(That's just a guess). Most of the time the sum produces identical
results, but for some numbers, it is different. The following
example, shows what happens:
set.seed(12233)
n <- 5
a <- runif(n,1,5)
e <- runif(n, 5*(1:n),10*(1:n))
s <- runif(1, 1.2, 4)
p <- runif(5, 3, 10)
x <- c(e[-5], (sum(e*p)-sum(e[-5]*p[-5]))/p[5])
u <- a^(1/s)*x^((s-1)/s)
dyn.load("sumF.so")
u[1] <- u[1]+.0001 ### If we do not add .0001, all differences are 0
s1 <- sum(u)
s2 <- .Fortran("sumf", as.double(u), as.integer(n), sf1=double(1),
sf2=double(1))[3:4]
s3 <- .C("sumc", as.double(u), as.integer(n), sC=double(1))[[3]]
s1-s2[[1]] ## R versus compiler sum() (Fortran)
[1] -7.105427e-15
s1-s2[[2]] ## R versus manual sum (Fortran
[1] -7.105427e-15
s1-s3 ## R Versus manual sum in C
[1] -7.105427e-15
s2[[2]]-s2[[1]] ## manual sum versus compiler sum() (Fortran)
[1] 0
s3-s2[[2]] ## Fortran versus C
[1] 0
My sumf and sumc are
subroutine sumf(x, n, sx1, sx2)
integer i, n
double precision x(n), sx1, sx2
sx1 = sum(x)
sx2 = 0.0d0
do i=1,n
sx2 = sx2+x(i)
end do
end
void sumc(double *x, int *n, double *sum)
{
int i;
double sum1 = 0.0;
for (i=0; i< *n; i++) {
sum1 += x[i];
}
*sum = sum1;
}
Can that be a bug? Thanks.
--
Pierre Chaussé
Assistant Professor
Department of Economics
University of Waterloo
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Discrepancy: R sum() VS C or Fortran sum
Hi all,
I found a discrepancy between the sum() in R and either a sum done in C
or Fortran for vector of just 5 elements. The difference is very small,
but this is a very small part of a much larger numerical problem in
which first and second derivatives are computed numerically. This is
part of a numerical method course I am teaching in which I want to
compare speeds of R versus Fortran (We solve a general equilibrium
problem all numerically, if you want to know). Because of this
discrepancy, the Jacobian and Hessian in R versus in Fortran are quite
different, which results in the Newton method producing a different
solution (for a given stopping rule). Since the solution computed in
Fortran is almost identical to the analytical solution, I suspect that
the sum in Fortran may be more accurate (That's just a guess). Most of
the time the sum produces identical results, but for some numbers, it is
different. The following example, shows what happens:
set.seed(12233)
n <- 5
a <- runif(n,1,5)
e <- runif(n, 5*(1:n),10*(1:n))
s <- runif(1, 1.2, 4)
p <- runif(5, 3, 10)
x <- c(e[-5], (sum(e*p)-sum(e[-5]*p[-5]))/p[5])
u <- a^(1/s)*x^((s-1)/s)
dyn.load("sumF.so")
u[1] <- u[1]+.0001 ### If we do not add .0001, all differences are 0
s1 <- sum(u)
s2 <- .Fortran("sumf", as.double(u), as.integer(n), sf1=double(1),
sf2=double(1))[3:4]
s3 <- .C("sumc", as.double(u), as.integer(n), sC=double(1))[[3]]
s1-s2[[1]] ## R versus compiler sum() (Fortran)
[1] -7.105427e-15
s1-s2[[2]] ## R versus manual sum (Fortran
[1] -7.105427e-15
s1-s3 ## R Versus manual sum in C
[1] -7.105427e-15
s2[[2]]-s2[[1]] ## manual sum versus compiler sum() (Fortran)
[1] 0
s3-s2[[2]] ## Fortran versus C
[1] 0
My sumf and sumc are
subroutine sumf(x, n, sx1, sx2)
integer i, n
double precision x(n), sx1, sx2
sx1 = sum(x)
sx2 = 0.0d0
do i=1,n
sx2 = sx2+x(i)
end do
end
void sumc(double *x, int *n, double *sum)
{
int i;
double sum1 = 0.0;
for (i=0; i< *n; i++) {
sum1 += x[i];
}
*sum = sum1;
}
Can that be a bug? Thanks.
--
Pierre Chaussé
Assistant Professor
Department of Economics
University of Waterloo
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Passing R code from webpage
If R and ghostcript are installed on the server that hosts you webpage, it is easy. All you need is a minimum of php. I started with doR (which I think does not exist anymore) and modified it. Even better, some people offers a solution for you. Here is a GPL licenced solution. http://steve-chen.net/document/r/r_php Le 2013-02-16 18:48, Matevz Pavlic a écrit : Hi all, Is there a way to pass R code from web page (html file) to do some statistics and than plot the output in web browser. I am looking forever at this, and cant find a way. Regards,m -- View this message in context: http://r.789695.n4.nabble.com/Passing-R-code-from-webpage-tp4658800.html Sent from the R devel mailing list archive at Nabble.com. __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] GPU Computing
Hi all, I am looking for a function similar to mclapply() that would work with GPU cores. I have looked at all possible packages related to GPU computing but they are mainly providing functionality for big dataset or big matrices. I use mainly mclapply to speed up simulations by running several simulations in parallel, which works nicely. Is it possible to do the same with a multicore GPU? I am planning to buy a tesla 2075, and I want to know before if it is something we can do with R. May be by modifying mclapply(). Thanks for your suggestions. -- *Pierre Chaussé* Assistant Professor Department of Economics University of Waterloo [[alternative HTML version deleted]] __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] kernapply.ts
I have a suggestion for kernapply for ts objects. When we choose the
option circular=F, the returned series don't have the correct dates. The
removed dates are all at the beginning instead of half at the beginning
and half at the end. It is particularly useful when we need to smooth
the series (or remove a trend using a filter) before estimating a model
(like in macroeconomics) or simply to plot the original series with the
smoothed one. Of course, there is always the option of doing it by hand
of the use circular=T and trim the series but I thought it would be
nicer that way.
Here is my suggestion (maybe not the nicest way to do it but it works)
kernapply.ts - function (x, k, circular = FALSE, ...)
{
if (!is.matrix(x))
{
y - kernapply.vector(as.vector(x), k, circular=circular)
ts (y, end=end(x), frequency=frequency(x))
}
else
y - apply(x, MARGIN=2L, FUN=kernapply, k, circular=circular)
if(circular)
ts (y, end=end(x), frequency=frequency(x))
else
{
y - as.ts(y)
t1 - tsp(x)[1]+(length(k[[1]])-1)/tsp(x)[3]
t2 - tsp(x)[2]-(length(k[[1]])-1)/tsp(x)[3]
tsp(y) - c(t1,t2,tsp(x)[3])
return(y)
}
}
--
*Pierre Chaussé*
Assistant Professor
Department of Economics
University of Waterloo
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] kernapply.ts
On 02/11/11 12:20 PM, Prof Brian Ripley wrote:
On Wed, 2 Nov 2011, Pierre Chausse wrote:
I have a suggestion for kernapply for ts objects. When we choose the
option circular=F, the returned series don't have the correct dates. The
That's a matter of opinion. A kernel is applied in the same way as an
MA filter, to historical data.
I understand for MA which is the weighted sum of past data but kernapply
does not average present and past values of the series X(t) but values
around it
ex. with length(k[[1]])=2
Smoothed(X_t) = k[[1]][2] X_{t-1} + k[[1]][1] X_{t} + k[[1]][2] X_{t+1}
which makes it natural to have the same date as X_t. Furthermore in the
kernapply.vector which is used for time series, the function returns the
following for circular=F
return (y[(1L+m):(n-m)])
In other words it removes the first ans last observations.
removed dates are all at the beginning instead of half at the beginning
and half at the end. It is particularly useful when we need to smooth
the series (or remove a trend using a filter) before estimating a model
(like in macroeconomics) or simply to plot the original series with the
smoothed one. Of course, there is always the option of doing it by hand
of the use circular=T and trim the series but I thought it would be
nicer that way.
Here is my suggestion (maybe not the nicest way to do it but it works)
kernapply.ts - function (x, k, circular = FALSE, ...)
{
if (!is.matrix(x))
{
y - kernapply.vector(as.vector(x), k, circular=circular)
ts (y, end=end(x), frequency=frequency(x))
}
else
y - apply(x, MARGIN=2L, FUN=kernapply, k, circular=circular)
if(circular)
ts (y, end=end(x), frequency=frequency(x))
else
{
y - as.ts(y)
t1 - tsp(x)[1]+(length(k[[1]])-1)/tsp(x)[3]
t2 - tsp(x)[2]-(length(k[[1]])-1)/tsp(x)[3]
tsp(y) - c(t1,t2,tsp(x)[3])
return(y)
}
}
--
*Pierre Chauss?*
Assistant Professor
Department of Economics
University of Waterloo
[[alternative HTML version deleted]]
--
*Pierre Chaussé*
Assistant Professor
Department of Economics
University of Waterloo
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] 0.5 != integrate(dnorm,0,20000) = 0
The warning about absolute error == 0 would not be sufficient
because if you do
integrate(dnorm, 0, 5000)
2.326323e-06 with absolute error 4.6e-06
We get reasonable absolute error and wrong answer. For very high upper
bound, it seems more stable to use Inf. In that case, another
.External is used which seems to be optimized for high or low bounds:
integrate(dnorm, 0,Inf)
0.5 with absolute error 4.7e-05
On 10-12-07 8:38 AM, John Nolan wrote:
I have wrestled with this problem before. I think correcting
the warning to absolute error ~= 0 is a good idea, and printing
a warning if subdivisions==1 is also helpful. Also, including
a simple example like the one that started this thread on the
help page for integrate might make the issue more clear to users.
But min.subdivisions is probably not. On the example with dnorm( ),
I doubt 3 subdivisions would work. The problem isn't that
we aren't sudividing enough, the problem is that the integrand
is 0 (in double precision) on most of the region and the
algorithm isn't designed to handle this. There is no way to
determine how many subdivisions are necessary to get a reasonable
answer without a detailed analysis of the integrand.
I've gotten useful results with integrands that are monotonic on
the tail with a self-triming integration routine
like the following:
right.trimmed.integrate- function( f, lower, upper, epsilon=1e-100,
min.width=1e-10, ... ) {
+ # trim the region of integration on the right until f(x) epsilon
+
+ a- lower; b- upper
+ while ( (b-amin.width) (f(b)epsilon) ) { b- (a+b)/2 }
+
+ return( integrate(f,a,b,...) ) }
right.trimmed.integrate( dnorm, 0, 2 ) # test
0.5 with absolute error 9.2e-05
This can be adapted to left trim or (left and right) trim,
abs(f(x)-c)epsilon,
etc. Setting the tolerances epsilon and min.width is an issue,
but an explicit discussion of these values could encourage people to
think about the problem in their specific case. And of course, none
of this guarantees a correct answer, especially if someone tries this
on non-monotonic integrands with complicated 0 sets. One could write
a somewhat more user-friendly version where the user has to specify
some property (or set of properties) of the integrand, e.g. right-tail
decreasing to 0, etc. and have the algorithm try to do smart
trimming based on this. But perhaps this getting too involved.
In the end, there is no general solution because any solution
depends on the specific nature of the integrand. Clearer messages,
warnings in suspicious cases like subdivisions==1, and a simple
example explaining what the issue is in the help page would help
some people.
John
...
John P. Nolan
Math/Stat Department
227 Gray Hall
American University
4400 Massachusetts Avenue, NW
Washington, DC 20016-8050
jpno...@american.edu
202.885.3140 voice
202.885.3155 fax
http://academic2.american.edu/~jpnolan
...
-r-devel-boun...@r-project.org wrote: -
To: r-devel@r-project.org, Prof Brian Ripleyrip...@stats.ox.ac.uk
From: Martin Maechler
Sent by: r-devel-boun...@r-project.org
Date: 12/07/2010 03:29AM
Subject: Re: [Rd] 0.5 != integrate(dnorm,0,2) = 0
Prof Brian Ripleyrip...@stats.ox.ac.uk
on Tue, 7 Dec 2010 07:41:16 + (GMT) writes:
On Mon, 6 Dec 2010, Spencer Graves wrote:
Hello:
The example integrate(dnorm,0,2) says it fails on many
systems.
I just got 0 from it, when I should have gotten either an error or
something
close to 0.5. I got this with R 2.12.0 under both Windows Vista_x64
and
Linux (Fedora 13); see the results from Windows below. I thought
you might
want to know.
Well, isn't that exactly what the help page says happens? That
example is part of a section entitled
## integrate can fail if misused
and is part of the illustration of
If the function is
approximately constant (in particular, zero) over nearly all its
range it is possible that the result and error estimate may be
seriously wrong.
yes, of course,
and the issue has been known for ``ages'' ..
..
..
but it seems that too many useRs are not reading the help
page carefully, but only browse it quickly.
I think we (R developers) have to live with this fact
and should consider adapting to it a bit more, particularly in
this case (see below)
Thanks for all your work in creating and maintaining R.
Best Wishes,
Spencer Graves
###
integrate(dnorm,0,2) ## fails on many systems
0 with absolute error 0
and this is particularly unsatisfactory for another reason:
absolute
[Rd] stats::kernel
Hi,
There is a small bug in the kernel() function. Everything is fine when
we use the format:
kernel(name,m,r)
but if we want the first argument to be a vector, which is useful is we
are interested in using a method not implemented in kernel(), the
default value of m is wrong. For example, if we do:
s - rep(1/11,6)
k - kernel(s)
we get the error message
Error in kernel(s) : 'coef' does not have the correct length
The problem is that the default value, which is not indicated in the
help file, violate the condition that length(coef) must be equal to
(m+1). Therefore, the first line of the function, which is:
function (coef, m = length(coef) + 1, r, name = unknown)
{
should be changed to
function (coef, m = length(coef) - 1, r, name = unknown)
{
bests
--
*Pierre Chaussé*
Assistant Professor
Department of Economics
University of Waterloo
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
