代写 R algorithm STAT 513/413: Lecture 3 The spirit of R

STAT 513/413: Lecture 3 The spirit of R
(Stand on the shoulders of experts)

We arrived to something – but it is not the end
1

R spirit?
Well, this aspect is not that easily encapsulated into few guidelines – we will rather strive all this course to get an idea what it is
But one thing we may start immediately with: avoid loops…
… think in terms of vectors/matrices, if possible
Another one, related to the previous one: use the code of experts
2

Vectorized
3

Finally, some comments perhaps
Seems like this is everything we can do about it
4

A tale of expert code I: floating point arithmetics
Floating-point arithmetics: numbers in R are represented as base ∗ 10exponent
which has inevitable consequences
> 0.000001*1000000
[1] 1
> x=0; for (k in (1:1000000)) x=x+0.000001
>x
[1] 1
> x=1000000; for (k in (1:1000000)) x=x+0.000001 >x
[1] 1000001
> x-1000000
[1] 1.000008
The moral here is: with floating-point arithmetics, adding works well if the added numbers are about of the same magnitude
5

A better algorithm thus does it
> x=0; for (k in (1:1000000)) x=x+0.000001; x=x+1000000 >x
[1] 1000001
> x-1000000
[1] 1
Yeah, but what to do in general? The solution seems to be: use addition programmed by experts
> sum
function (…, na.rm = FALSE) .Primitive(“sum”)
> sum(c(1000000,rep(0.000001,1000000)))
[1] 1000001
> sum(c(1000000,rep(0.000001,1000000)))-1000000
[1] 1
6

Vectorization alone does it not
> rep(1,1000001) %*% c(1000000,rep(0.000001,1000000))-1000000
[,1]
[1,] 1.000008
> crossprod(rep(1,1000001),c(1000000,rep(0.000001,1000000)))-1000000
[,1]
[1,] 1.000008
7

Appendix. More modi operandi in R: functions
• function: the input can be varied in a better way than was reediting the script, and the variables inside the function do not mess up in your working environment (scoping)
line <- function(n,a,b,u=1,v=10,s=1) ### plots n points approximately following a line with given ### intercept and slope, and normal error controled by sigma ### x-points sampled from uniform(lower,upper) { x=runif(n,u,v) y=a+b*x+rnorm(n) plot(x,y) } All this process enables you to vary input - first in script, then in function - and thus get some more confidence that the whole concoction does the right thing However, for this course we are just fine with scripts - although successful scripts can be easily upgraded to functions, and those are allowed as well 8