# 0 means skip, 1 means do, but don't print figs, 2 means
# also print figs to file
doprob1 = 1; #This never go written, so we renumbered
doprob2 = 1; # Goes with prob 1 in pset pdf file
doprob3 = 1; # Goes with prob 2 in pset pdf file
doprob4 = 1; # ...
doprob5 = 1; # ...
doprob6 = 1; #...

setpdf = function(fname,w,h)
{
    pdf(file=fname, width=w, height=h)
}
endpdf = function()
{
    dev.off()
}

if (doprob1)
{
    #NEVER WROTE THIS PROBLEM
    
    x = 1
    #rows H_R H_Y H_G columns red sensor, yellow sensor,
    #green sensor
    #probtable=    [
    #    s1 s2 s3     000 001 010 011 100  101  110  111
    #    Hr     1   .6  .8;  0    0   0   0  .08  .32  .12  .48
    #    Hy     .7  .8  .8;  
    #    Hg     .5  .6  .9]
    
    #s2
    #s3
    #prior = [1/3 1/3 1/3]
    #
    #[0 0 0]
    #0
    #Hr
    #Hy
    #Hg
}

if (doprob2)
{
    lw = 2;
    a = c(1, 10, 50, 500, 30)
    b = c(1,10, 50, 500, 70)
    maxind = 4 #index of maximum indices
    col = c(rgb(0,1,0), rgb(0,1,1), rgb(1,1,0), rgb(1,0,0), rgb(0,0,1))
    nprior = 5

    print('2a')
    lbf = 250*log(.5) - 140*log(140/250) - 110*log(110/250)
    bf = exp(lbf)
  
    print('2b')
    if (doprob2 == 2)
    {      
      setpdf('../pdftex/img/ps6-2b.pdf', 5,4)
      par(mar=c(4,2,1,1))
    }
    th = seq(0,1,.0005)
    n = 0;
    m = 0;
    xlab = expression(theta)
    j = maxind
    plot(th,dbeta(th,a[j]+n,b[j]+m), type='l',col=col[j], lwd=lw, 
         xlab=xlab,ylab='',cex.lab=2)
    for (j in 1:nprior)
    {
        if (j == maxind)
          next
        lines(th,dbeta(th,a[j]+n,b[j]+m), col=col[j], lwd=lw)
    }

    abline(v=.5,col='red',lty='dashed')
    #m = {'Beta(1,1)', 'Beta(10,10)', 'Beta(50,50)', ...
    legend(.55,25, 
           c('Beta(1,1)', 'Beta(10,10)', 'Beta(50,50', 'Beta(500,500)', 'Beta(30,70)'),
             lty=c(1,1), lwd=2, col=col)
    if (doprob2 == 2)
      endpdf()

  print('2c')
  n = 0;
  m = 0;
  for (j in 1:nprior)
    print(1- pbeta(.5,a[j]+n,b[j]+m))
  
# Problem 2d ----
  print('2d')
  if (doprob2 == 2)
  {
      setpdf('../pdftex/img/ps6-2d.pdf',5,4)
      par(mar=c(4,2,1,1))
  }
  th = seq(0,1,.001)
  n = 140;
  m = 110;
  j = maxind
  r = a[j]+n;
  s = b[j]+m;    
  xlab = expression(theta)
  plot(th,dbeta(th,r,s), type='n', col=col[j],lwd=lw,xlab=xlab,ylab='',
       cex.lab = 2)
  for (j in 1:5)
  {
      r = a[j]+n;
      s = b[j]+m;
      lines(th,dbeta(th,r,s), col=col[j],lwd=lw)
      v = dbeta(.6,r,s);
      st = sprintf('Beta(%d,%d) at .6 = %f', r, s, v);
      print(st)
  }
  abline(v=140/250, col='red', lty='dashed')
  abline(v=.5, col='red', lty='dashed')
  legend(.6,25, 
         c('Beta(1,1)', 'Beta(10,10)', 'Beta(50,50', 'Beta(500,500)', 'Beta(30,70)'),
         lty=c(1,1), lwd=2, col=col)
  if (doprob2 == 2)
    endpdf()

  print('2e')
  n = 140;
  m = 110;
  for (j in 1:nprior)
  {
    r = a[j]+n;
    s = b[j]+m;
    v = 1- pbeta(.5,r,s)
    st = sprintf('1 - pbeta(0.5,%d,%d)  = %f', r, s, v);
    print(st)
  }

  print('2f')
  n = 140;
  m = 110;
  for (j in 1:nprior)
  {
    r = a[j]+n;
    s = b[j]+m;
    v = dbeta(.56,r,s)/dbeta(.5,r,s)
    st = sprintf('dbeta(.56,%d,%d)/dbeta(.5,%d,%d)  = %f', r, s, r,s,v);
    print(st)
  }
dbeta(.56,a[j]+n,b[j]+m)/dbeta(.5,a[j]+n,b[j]+m)
}

if (doprob3)
{
    d = c(12, 10, 11, 4, 11)
    al = 10;
    bl = 15;
  
    loglikalice = sum(d)*log(al) - length(d)*al; #ignore factorials
    loglikbob = sum(d)*log(bl) - length(d)*bl;   #ignore factorials
    
    bf = exp(loglikalice - loglikbob)
    odds = bf*.1
    print(odds)
}

if (doprob4)
{
    if (doprob4 == 2)
        setpdf('../pdftex/img/ps6-4a.pdf')
    lw = 2;
    print('4a: figure 1')
    x = .2;
    c1 = -1/log(x);
    col = 'red';
    xlab = expression(theta)
    plot(c(0,1), c(0,3.5), type='n', xlab=xlab, ylab='')
    lines(c(x,x),c(0,c1/x), col=col, lwd=lw)
    th = seq(x,1,.01)
    lines(th, c1/th, col=col, lwd=lw)
    s = sprintf('posterior(.3) for x=.2: %f', c1/.3);
    print(s)
    
    x = .5;
    c1 = -1/log(x);
    col = 'blue';
    lines(c(x,x),c(0,c1/x), col=col, lwd=lw)
    th = seq(x,1,.01)
    lines(th, c1/th, col=col, lwd=lw)

    s1 = expression(paste('f(',theta, '|x = .2'))
    s2 = expression(paste('f(',theta, '|x = .5'))
    legend(.6,3.5,c(s1,s2),col=c('red','blue'), lty=1,lwd=3)
    
    if (doprob4 == 2)
        endpdf()

    print('4b: figure 2')
    if (doprob4 == 2)
    {
      setpdf('../pdftex/img/ps6-4b.pdf',5,4)
      par(mar=c(4,2,1,1))
    }
    
    xm = .5;
    n = 2;
    c1 = (n-1)/(xm^(1-n)-1);
    col = 'red';
    plot(c(0,1),c(0,10), col=col, lwd=lw, type='n', 
         xlab=expression(theta), ylab='')
    
    #Draw legend first so it doesn't cover the curves
    s1 = 'data = (.1, .5) or (.5,.5)'
    s2 = 'data = (.1, .2, ..., .5)'
    legend(.0,10,c(s1,s2),col=c('red','blue'), lty=1,lwd=3)
    
    th = seq(xm,1,0.01)
    lines(th, c1/(th^n), col=col, lwd=lw)
    s = sprintf('posterior(.7) for data (.1,.5): %f', c1/(.7)^n);
    print(s)

    xm = .5;
    n = 5;
    c1 = (n-1)/(xm^(1-n)-1);
    col = 'blue';
    lines(c(0,xm,xm),c(0,0,c1/xm^n), col=col, lwd=lw)
    th = seq(xm,1,.01)
    lines(th, c1/(th^n), col=col, lwd=lw)
    s = sprintf('posterior(.7) for data (.1,.2,.3,.4,.5): %f', c1/(.7)^n);
    print(s)

    if (doprob4 == 2)
        endpdf()
  
    print('4c')
    xm = .5;
    n = 5;
    c1 = (n-1)/(xm^(1-n)-1);
    predprob = (-.5*c1/5)*(1-.5^(-5))  

    print('4d: figure 3')
    if (doprob4 == 2)
    {
      setpdf('../pdftex/img/ps6-4d.pdf',5,4)
      par(mar=c(4,2,1,1))
    }
    
    col = 'red';
    plot(c(0,1),c(0,1), type='l', col=col, lwd=lw,xlab='x', ylab='')
    
    col = 'blue';
    x = seq(0,.01,.0001)
    lines(x,-(1-x)/log(x), col=col, lwd=lw)
    x = seq(.01,.998,.002)
    lines(x,-(1-x)/log(x), col=col, lwd=lw)

    s1 = 'MAP'
    s2 = 'Conditional expectation'
    legend(0.3,0.25,c(s1,s2),col=c('red','blue'), lty=1,lwd=3)
  
    if (doprob4 == 2)
        endpdf()
}

if (doprob5)
{
    print('5a')
    muprior = 5;
    varprior = 9;
    vartheta = 4;
    n=1;
    xbar = 6;
    a = 1/varprior;
    b = n/vartheta;
    mupost = (a*muprior + b*xbar)/(a+b)
    varpost = 1/(a+b)

    print('5b: figure 1')
    if (doprob5 == 2)
    {
      setpdf('../pdftex/img/ps6-5b.pdf',5,4)
      par(mar=c(4,2,1,1))
    }
    
    lw = 2;
    muprior = 5;
    varprior = 9;
    vartheta = 4;
    n=4;
    xbar = 6;
    a = 1/varprior;
    b = n/vartheta;
    mupost = (a*muprior + b*xbar)/(a+b)
    varpost = 1/(a+b)

    th = seq(-3,13,.01)
    col = 'red';
    plot(c(-5,15),c(0,.5), type='n',xlab=expression(theta), ylab='')
    
    #Draw legend first to avoid covering curves
    s1 = 'Prior N(5,9)'
    s2 = 'Posterior N(5.9,0.9)'
    legend(-5.5,0.5,c(s1,s2),col=c('red','blue'), lty=1,lwd=3)
    
    lines(th, dnorm(th, muprior,sqrt(varprior)), col=col, lwd=lw)
        
    col = 'blue';
    lines(th, dnorm(th, mupost,sqrt(varpost)), col=col, lwd=lw)
    
    if (doprob5 == 2)
        endpdf()

    print('5d')
    muprior = 100;
    varprior = 225;
    vartheta = 100;
    n=1;
    a = 1/varprior;
    b = n/vartheta;
    xbar = 80;
    mupost = (a*muprior + b*xbar)/(a+b)
    varpost = 1/(a+b)
    
    xbar = 150;
    mupost = (a*muprior + b*xbar)/(a+b)
    varpost = 1/(a+b)
}

if (doprob6)
{
    data = c(1, 0, 1, 1, 1)
    bfactors = c((3/4)/(5/6), (1/4)/(1/6))
    priorodds = 1
    odds = priorodds;
    for (dataval in data)
    {
        ind = dataval + 1;  #dataval is 0 or 1, column index is 1 or 2
        odds = odds*bfactors[ind]
    }
}