0. Bonev's info.
1. Projects.
2. Bias in sample variance?
3. Making R packages.
4. Calling R from C.

0. Bonev's Python and object oriented programming info.
Boyan Bonev's email is bonev@ucla.edu.
His slides are at 
http://www.stat.ucla.edu/~boyan/python_stats202a.zip
http://www.stat.ucla.edu/~boyan/pythonpdfs_stats202a.zip
and http://www.slideshare.net/nowells/introduction-to-python-5182313

1. Projects. 

See lecture 10 for my rules. 
Rule 4 is below. 

Rule 4: Send me a pdf or powerpoint version of your slides by 6pm the night 
before your talk, to frederic@stat.ucla.edu. That way, I can set up the 
talks in order ahead of time and we won't have to 
waste time in class waiting for each person to connect their laptop to the 
projector. About 7 or 8 slides seems right, though it's fine with me if 
you want fewer or more.

Do not feel the need to show us ALL your results.

2. Bias in sample variance?

The sample variance s^2 = · (Xi - Xbar)^2 / (n-1).
s^2 is an unbiased estimate of sigma^2, regardless of n and the 
distribution of Xi. 
Last time we investigated the bias in s as an estimate of sigma,
supposing the random variables Xi are uniform (0,1).
What would it look like if we repeated this, 
but looked at s^2 instead of s?

In sd3.c,

#include <R.h>
#include <Rmath.h>

double s2(double *x, int j) {
    int i;
    double a1, xbar;
    
    xbar = 0.0;
    for(i = 0; i < j; i++) {
	xbar += x[i];
    }
    xbar = xbar / j;
    a1 = 0.0;
    for(i = 0; i < j; i++) {
	a1 += pow(x[i] - xbar,2);
    }
    a1 = a1/(j-1); // I changed this from a1 = pow(a1/(j-1), 0.5).
    return(a1);
}

double mean2(double *x, int j) {
    int i;
    double xbar;
    
    xbar = 0.0;
    for(i = 0; i < j; i++) {
	xbar += x[i];
    }
    xbar = xbar / j;
    return(xbar);
}

void bias4 (int *nmax, int *m, double *b){
    /* The output will be b, the bias in s2 for each n,
         from n=2 up to nmax. */
    int i,j,k;
    double x[*nmax];
    double s[*m];
    
    GetRNGstate();
    for(i = 2; i < *nmax+1; i++){   
	for(j = 0; j < *m; j++){
	    for(k = 0; k < i; k++){
		x[k] = runif(0,1);
	    }
	    s[j] = s2(x, i); // I changed this from sd2 to s2.
	}
	b[i-1] = mean2(s, *m) - 1.0/12.0; // I changed pow(12,-.5) to 1.0/12.0.
	if(i/1000 == floor(i/1000)) Rprintf("%i ",i);
    }
    PutRNGstate();
}

In R, 

  system("R CMD SHLIB sd3.c")
  dyn.load("sd3.so")
  nmax = 100
  m = 100000
  a = .C("bias4",as.integer(nmax), as.integer(m),b = double(nmax))
  bias1 = a$b
  plot(c(2:nmax),bias1[2:nmax],xlab="n",ylab="bias", ylim=c(-.003,.001),
         main="bias in SD of n uniforms", type="l")
  abline(h=0,lty=2)

3. Making R packages. 
See buildingRPackages.pdf.

library(holdem)
name1 = c("gravity","tommy","ursula","timemachine","vera","william","xena")
decision1 = list(gravity, tommy, ursula, timemachine, vera, william, xena) 
tourn1(name1, decision1, myfast1 = 0, t1=0, t2=0.5)
tourn1(name1, decision1, myfast1 = 2)

n = choose(44,3); result = rep(0,n); a1 = c(8,22,34,35,48,49,50,51)
a2 = c(1:52)[-a1]; i= 0
for(i1 in 1:42){for(i2 in ((i1+1):43)){for(i3 in ((i2+1):44)){
    flop1 = c(a2[i1],a2[i2],a2[i3]); flop2 = switch2(flop1)
    b1 = handeval(c(10,10,flop2$num),c(3,2,flop2$st))
    b2 = handeval(c(9,9,flop2$num),c(3,1,flop2$st))
    b3 = handeval(c(12,10,flop2$num),c(4,4,flop2$st))
    b4 = handeval(c(13,11,flop2$num),c(4,4,flop2$st))
    i = i+1;  if(b1 > max(b2,b3,b4)) result[i] = 1}}
    cat(i1)}
sum(result > 0.5)

n = 10000; a = rep(0,n)
for(i in 1:n){
    x1 = deal1(2)
    b1 = handeval(c(x1$plnum1[1,],x1$brdnum1), c(x1$plsuit1[1,],x1$brdsuit1))
    b2 = handeval(c(x1$plnum1[2,],x1$brdnum1), c(x1$plsuit1[2,],x1$brdsuit1))
    if(min(b1,b2) > 4000000) a[i] = 1
    if(i/1000 == floor(i/1000)) cat(i)
}
sum(a>.5)

4. Calling R from C.

Notes on this come from Phil Spector's lecture notes at UC Berkeley.

The composite version of Simpson's rule is a rough way to approximate 
an integral of some function, f.
It goes back to English mathematician Thomas Simpson (1710Ð1761), who 
approximated the integral from a to b of f(x) dx by 
[(b-a)/6] [f(a) + 4f((a+b)/2) + f(b)].
This is the approximation you'd get by just dividing the interval (a,b) 
into 6 pieces, and using the left endpoint for the first piece, the right 
endpoint for the last piece, and the middle, (a+b)/2, for the other 4 
pieces.

For the composite Simpson's rule, you divide the interval not into 
6 pieces but n pieces, where n is even, and get
the integral from a to b of f(x) dx ~
  h/3 [f(a) + 4f(a+h) + 2f(a+2h) + 4f(a+3h) + 2f(a+4h) + ... + 4f(b-h) + f(b)],
where h = (b-a)/n.
Note that the coefficients 4 and 2 alternate.
The error in Simpson's rule is bounded in absolute value by 
h^4 (b-a) max[f^4(x); x in (a,b)]/180.

Suppose we want to write this function f in R, but run the composite Simpson's
rule in C, and call the C function from R.
The key component is a C function called call_R. call_R takes as one of its 
arguments a pointer to a function, which is passed into the C environment 
as a list. The number and type of the arguments to f are also passed to call_R. 
As usual, when calling C from R, the results of the C function in the end are 
passed back to R not as a returned value, but through the argument list. 

It is typically easiest to divide the C program into three parts.
1) simp is the algorithm to implement the composite Simpson's rule, 
but it doesn't deal with issues pertaining to the interface between C and R.
2) sfunc is just a shell that uses call_R to have R evaluate the function f.
3) dosimp is just the interface, which is called by R and which calls simp.

In simp.c,

  static char *sfunction; 

 double simp(double(*func)(),double start,double stop,long n)
{
    double mult,x,t,t1,inc; 
    long i; 
    inc = (stop - start) / (double)n; 
    x = start;
    t = func(x); 
    mult = 4.0; 
    for(i=1; i<n; i++){
	x += inc; 
	t += mult * func(x); 
	if(mult == 4.0) mult = 2.0;
	else mult = 4.0;
    }
    t += func(stop); 
    return(t * inc / 3.0); 
}

  double sfunc(double x) {
    char *modes[1]; 
    char *arguments[1]; 
    double *result; 
    long lengths[2]; 
    lengths[0] = (long)1; 
    arguments[0] = (char*)&x; 
    modes[0] = "double"; 
    call_R(sfunction,(long)1,arguments,modes,lengths,
	   (char*)0,(long)1,arguments); 
	   result = (double*)arguments[0]; 
	   return(*result); 
  }

  dosimp(char** funclist,double *start,double *stop,long *n,double *answer) {
    double sfunc(double); 
    double simp(double(*)(),double,double,long); 
    sfunction = funclist[0]; 
    *answer = simp((double(*)())sfunc,*start,*stop,*n); 
  } 
	   
The function simp implements the composite Simpson's rule.

The function dosimp is what R calls.
dosimp simply extracts the function pointer from the list passed into the C 
environment, and loads it into a static pointer to make it available to the 
function sfunc, 
and then calls the function simp which actually does the integration. 

sfunc uses the call_R function to pass arguments to, and receive results from, R. 
call_R's first argument is the pointer to the function which was passed to dosimp from R. 
The second argument is the number of arguments to be passed to the R 
function, which in this case is 1. (long 1 just means the integer 1.)
Next is an array of pointers to the arguments themselves, each cast as char*. 
Note that this is an array of pointers, i.e. a pointer to pointers.
Similarly, the modes argument is an array of pointers to character strings 
that tell the R function how to treat each argument, i.e. as doubles, 
integers, etc. The most common modes are double and character. 
As in all functions external to R, the arguments themselves are pointers 
to the values they represent, not the values themselves. 
The next argument is an array of integers (longs) giving the length of 
each of the arguments passed; in this case, lengths[0] has been set to 1. 
The next argument, which is cast as a null pointer in this example, 
is a pointer to an array of names of the arguments being passed 
to the R function; usually an explicit pointer to names will not 
be required. The next argument is the number of results which the 
R function will return; by returning a list, an R function can return 
more than one value into the C environment. 
The final argument is an array of pointers to where the results will be returned. 

Now, in R,
  system("R CMD SHLIB simp.c")
  dyn.load("simp.so")
  mysimp = function(func, start, stop, n = 10) {
    if(trunc(n/2) - n/2 != 0) stop("Number of points must be even.")
    t=0 
    .C("dosimp",list(func),as.double(start),
    as.double(stop),as.integer(n), ans = as.double(t))$ans
  }
  f1 = function(x){x^4 - 2*x^2}
  mysimp(f1,1,4,10)  ## compare with 4^5/5 - 1/5 - 2*4^3/3 + 2/3 = 162.6.
  mysimp(f1,1,4,100)

Note the use of list(func). This passes the function pointer to the C routine.