Solution of Homework Set 5, STAT 705,   Fall 2017.
=====================================================

## a(i)-(iii)

gfunc = function(x,y,c)  (3*x^2 + c*x*y + 5*y^3)/(1+x^2 + y^2)^2

MonteC.int = function(a,b,c, N=1000)  {
      xval = runif(N,-1,1)*a
      yval = runif(N,-1,1)*b
      gval = gfunc(xval,yval,c)
      list(intg = 4*a*b*mean(gval), SE = 4*a*b*sd(gval)/sqrt(N)) }

gint.innr = function(xvec,b,c, reltol=1e-6) {
    ### xvec is vector of arguments for \int_b^b g(x,y) dy 
    M = length(xvec)
    out = numeric(M)
    for(i in 1:M) out[i] = integrate(function(y) 
           gfunc(xvec[i],y,c),-b,b, rel.tol=reltol)$val
    out }

num.int = function(a,b,c,reltol=1e-6) 
    integrate(gint.innr,-a,a,rel.tol=reltol, b=b, c=c)[c(1,2)]


> MonteC.int(1,2,3,1e6)
$intg
[1] 1.529475         ### CI 1.529475 + c(-1,1)*0.0082705
                     ###   =   (1.521204, 1.537745)
$SE
[1] 0.008270531

> num.int(1,2,3,1e-5)
$value
[1] 1.531232

$abs.error
[1] 2.029205e-06

> num.int(1,2,3,1e-6)    ## same as with reltol=1e-9
$value
[1] 1.531232

$abs.error
[1] 1.700009e-14   ## not a real estimate of error because does not 
                  ## include errors of inner integral
                  ### but still seems to be very accurate

>  MonteC.int(1,2,3,3e7)
$intg
[1] 1.532385

$SE
[1] 0.001510567


### (b) 

lgsreg = function(XY, tol=1e-8) {
    ### first code the log-likelihood as function of a,b
    xvec = XY[,1]; yvec = XY[,2]
    nloglik = function(ab) -sum(yvec*(ab[1]+ab[2]*xvec)-
          log(1+exp(ab[1]+ab[2]*xvec)))
    tmp = optim(c(0,0), nloglik, method="BFGS", 
          hessian=T, control=list(reltol=tol))
    list(coef=tmp$par, SE = sqrt(diag(solve(tmp$hess)))) }

> xtmp = runif(100)
> ytmp = rbinom( 100, 1, plogis(-2 + 1.5*xtmp) )
> sum(ytmp)
[1] 25

> lgsreg(xtmp,ytmp)
$coef
[1] -2.246961  2.115283

$SE
[1] 0.5448664 0.8355406

> tmp = glm(ytmp ~ xtmp, family=binomial)
> summary(tmp)$coef[,1:2]
             Estimate Std. Error
(Intercept) -2.245969  0.5447548
xtmp         2.113709  0.8354336

> lgsreg(xtmp,ytmp,1e-10)
$coef
[1] -2.245969  2.113709

$SE
[1] 0.5447548 0.8354337

> abmat = array(0, c(1000,2))
  for(i in 1:1000) {
      Xv = runif(100)
      abmat[i,] = lgsreg(cbind(Xv, rbinom(100,1,
           plogis(-1 + 2*Xv))))[[1]]  
  }

> hist(abmat[,2], nclass=20, prob=T)
  curve(dnorm(x,mean(abmat[,2]), sd(abmat[,2])),
          add=T, col="blue")

### This last line not needed for HW solution