1. Projects.
2. Bias in sample variance?
3. Making R packages.

1. Projects. 

See lecture 11 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.

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.

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)