Solution to Homework 18
=======================

lowbd =  c(-3,-2,0,-1,-10,3,-8,-5,-2,-4,-1,-4,0,2,-5,-2)

### carr   16x7 matrix as given on web-page

X0 = rep(1,7)

Xtrans = function(Xini, carr, lowbd) {
   newbd = lowbd-c(carr %*% Xini)
   newv = rnorm(length(Xini))
   newv = newv/sqrt(sum(newv^2))
   dotp = c(carr %*% newv)
   indneg = (dotp < 0)
   tmax = min(newbd[indneg]/dotp[indneg])
   indpos = (dotp > 0)
   tmin = max(newbd[indpos]/dotp[indpos])
   list(Xout=Xini + runif(1,tmin,tmax)*newv, tmin=tmin, 
          tmax=tmax, newv=newv) }

tmpout=Xtrans(X0,carr,lowbd)
> lowbd-c(carr %*% tmpout$Xout)
 [1] -1.18998989 -1.72226423 -1.32168980 -1.11512491 -2.40024827 -0.04270208
 [7] -1.81185923 -1.06940398 -0.86528965 -0.84559418 -0.69235185 -2.64084950
[13] -1.98084622 -0.38746452 -2.15214309 -1.58421957

unlist(tmpout[2:3])
      tmin       tmax
-0.4230027  0.3601812

### Now we will create a sequence of 1000 iterates, after
###    pre-iterating 200

X1=tmpout$Xout
for(i in 1:200) X1 = Xtrans(X1,carr,lowbd)$Xout
Xmat = array(0, c(1000,7))
for(i in 1:1000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }

### To check that all of these are feasible points:

> summary( c(outer(lowbd,rep(1,1000)) - carr %*% t(Xmat) ) )
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-7.4190000 -2.5880000 -1.3900000 -1.8150000 -0.6853000 -0.0002688

### Feasible meant that each of the values in the 16 x 1000 array
    created was negative, so OK.

### Now need some checks of convergence.

#### >>>>> #### (First Check)
w = runif(7)
w = w/sqrt(sum(w*w))
[1] 0.4573 0.4842 0.1007 0.0294 0.4725 0.4072 0.3955

>  range(c(Xmat %*% w))
[1] 1.392719 3.206320

vec1=hist(c(Xmat %*% w), breaks=seq(0.85,3.45, by=0.1), prob=T)$counts
> vec1
 [1]   0   0   0   0   0   3  17  33  86  97  80  86 145 101  73  76  59  51  39
[20]  22  19   7   5   1   0   0

### Now pre-iterate for 5000, and try again.

X1 = Xmat[1000,]
for(i in 1:5000) X1= Xtrans(X1,carr,lowbd)$Xout
for(i in 1:1000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }
> range(c(Xmat %*% w))
[1] 1.460378 3.093886        ### Narrower range

vec2=hist(c(Xmat %*% w), breaks=seq(0.85,3.45, by=0.1), prob=T)$counts
 > vec2
 [1]   0   0   0   0   0   0   5  38  48  57  38  78  75  85  87 108  88  74  75
[20]  67  43  29   5   0   0   0

### Different shape !! But it is hard to know what the right shape is:
### let's create  a matrix and then print 10 different simulated
### ouptput counts for blocks of 1000, separated by `pre-iterated'
### blocks of 5000. We print the middle 14 categories corresponding
### to values 1.45--2.85.

CountMat = array(0,c(10,26))
RangMat = array(0,c(10,2))
CountMat[1,]=vec1; CountMat[2,]=vec2; 
RangMat[1,] =c(1.293272, 3.120429)
RangMat[2,] =c(1.276512, 2.739910)

for(j in 1:8) {
  X1 = Xmat[1000,]
  for(i in 1:5000) X1= Xtrans(X1,carr,lowbd)$Xout
  for(i in 1:1000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }
  proj = c(Xmat %*% w)
  RangMat[j+2,] = range(proj)
  CountMat[j+2,] = hist(proj, breaks=seq(0.85,3.45, by=0.1), prob=T)$counts }

round(RangMat,4)
        [,1]   [,2]
 [1,] 1.2933 3.1204
 [2,] 1.2765 2.7399
 [3,] 1.6365 3.1555
 [4,] 0.8659 3.1328
 [5,] 1.5188 3.2173
 [6,] 1.2595 3.1001
 [7,] 1.4764 3.0157
 [8,] 1.5336 3.1443
 [9,] 1.2628 3.1480
[10,] 1.4747 3.2455       ### Nothing too obvious except that the
 ###  variation among the first 5 is much more than among the second 5

> c(var(RangMat[1:5,1]),var(RangMat[6:10,1]))
[1] 0.08713050 0.01695485
> c(var(RangMat[1:5,2]),var(RangMat[6:10,2]))
[1] 0.036102164 0.006957089

CountMat[,7:20]

  17  33   86   97   80   86  145  101   73  76  59  51  39  22
  5   38   48   57   38   78   75   85   87 108  88  74  75  67
  0    2   28   18   43   27   58   83  120  99  91  96  85  97
 46   38   19   38   34   51   48   69  104 126 132  68  53  27
  5    6   15   19   51   62   78  129  105  77  86  97 103  82
 21   29   29   45   61   88   96   84   54  96  89  69  94  62
 13   22   49  101   86  116   92   84  107  59  86  98  27  40
  3   15   29   21   34   41   50   79   99 118 146 125 103  65
 67   46   37   75  128  122   88   71   59  56  59  55  44  24
 16   11   16   25   37   84   81   70  103 120 142 100  67  55
     

### Does not tell so much except that there is a lot of variability!

#### >>>>> #### (Second Check)
### Now re-do the whole thing, but overlay density estimates, with a
###  different w  for  6  different  blocks of 2000  separated by 3000

w = runif(7); w = w/sqrt(sum(w*w))
> round(w,4)
[1] 0.2287 0.4977 0.4973 0.3396 0.5524 0.1187 0.1346
> range(Xmat%*%w)
[1] 1.710972 3.463750

plot(c(Xmat%*%w), type="n", xlab="Projection on w", ylab="Density",
   xlim=c(1.7,3.5), ylim=c(0,1.2), main=paste("Overlaid Density", 
   " Estimates \n 6 Blocks of 2000 sep by 3000",sep=""))

Xmat = array(1, c(2000,7))
for(j in 1:6) {
   X1 = Xmat[2000,]
   for(i in 1:3000) X1= Xtrans(X1,carr,lowbd)$Xout
   for(i in 1:2000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }
   proj=c(Xmat %*% w)
   tmpdns = density(proj)
   lines(tmpdns$x,tmpdns$y, lty=j) }
legend(locator(), legend=paste("Blk",1:6), lty=1:6)
### Looks OK, cannot easily see glaring differences 

#### >>>>> #### (Third Check) not done here

Overlay  acf or spectral (method="ar") plots for separate blocks of
MCMC output converted into scalar series, to see that the blocks are
distributed more or less trhe same ,just as we did above for 
densities.


#### >>>>> #### (Fourth Check, confirming uniformity on a box)

## Find a box around rep(1,7) contained in the region. If it is not
   too small, then:

### Divide this box into 2^7 = 128 pieces by conditions that each
    coord is <1 or > 1, using "sign" function

### Now for each of a very large number of MCMC generated vectors
####   (1) check that all coord's within .074 of 1, discard if not.
####   (2) for those within range, create vector of signs,then takes
###        2*s-1 is each coord and inner product with 2^(0:6) and add
###        1. This gives numerical category 1:128 which should all be
###         equally likely.

#### This calculation will do two things:
####      (i)  check uniformity within the box
####      (ii) if we compare the total number in the box compared to
####           the total number N simulated, we get an estimate of 
####           the volume of the region, whose reciprocal is the
####           normalizing constant for the density !

Counts = numeric(128)
a0 = 2^(0:6)
## assume the box half-lengths (not nec equal) are given in vector "box"

## We start with the last MCMC iterate X1:  then

for(i in 1:1e5) {
        X1= Xtrans(X1,carr,lowbd)$Xout
        if(max(abs(X1-1))<box) {
           m = sum((sign(X1-1)+1)*a0)/2 + 1
           Counts[m] = Counts[m]+1 } } 

### For additional details, see the related text-file 
###    HW18slnaux.txt  in this Directory.Solution to Homework 18
=======================

lowbd =  c(-3,-2,0,-1,-10,3,-8,-5,-2,-4,-1,-4,0,2,-5,-2)

### carr   16x7 matrix as given on web-page

X0 = rep(1,7)

Xtrans = function(Xini, carr, lowbd) {
   newbd = lowbd-c(carr %*% Xini)
   newv = rnorm(length(Xini))
   newv = newv/sqrt(sum(newv^2))
   dotp = c(carr %*% newv)
   indneg = (dotp < 0)
   tmax = min(newbd[indneg]/dotp[indneg])
   indpos = (dotp > 0)
   tmin = max(newbd[indpos]/dotp[indpos])
   list(Xout=Xini + runif(1,tmin,tmax)*newv, tmin=tmin, 
          tmax=tmax, newv=newv) }

tmpout=Xtrans(X0,carr,lowbd)
> lowbd-c(carr %*% tmpout$Xout)
 [1] -1.18998989 -1.72226423 -1.32168980 -1.11512491 -2.40024827 -0.04270208
 [7] -1.81185923 -1.06940398 -0.86528965 -0.84559418 -0.69235185 -2.64084950
[13] -1.98084622 -0.38746452 -2.15214309 -1.58421957

unlist(tmpout[2:3])
      tmin       tmax
-0.4230027  0.3601812

### Now we will create a sequence of 1000 iterates, after
###    pre-iterating 200

X1=tmpout$Xout
for(i in 1:200) X1 = Xtrans(X1,carr,lowbd)$Xout
Xmat = array(0, c(1000,7))
for(i in 1:1000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }

### To check that all of these are feasible points:

> summary( c(outer(lowbd,rep(1,1000)) - carr %*% t(Xmat) ) )
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-7.4190000 -2.5880000 -1.3900000 -1.8150000 -0.6853000 -0.0002688

### Feasible meant that each of the values in the 16 x 1000 array
    created was negative, so OK.

### Now need some checks of convergence.

#### >>>>> #### (First Check)
w = runif(7)
w = w/sqrt(sum(w*w))
[1] 0.4573 0.4842 0.1007 0.0294 0.4725 0.4072 0.3955

>  range(c(Xmat %*% w))
[1] 1.392719 3.206320

vec1=hist(c(Xmat %*% w), breaks=seq(0.85,3.45, by=0.1), prob=T)$counts
> vec1
 [1]   0   0   0   0   0   3  17  33  86  97  80  86 145 101  73  76  59  51  39
[20]  22  19   7   5   1   0   0

### Now pre-iterate for 5000, and try again.

X1 = Xmat[1000,]
for(i in 1:5000) X1= Xtrans(X1,carr,lowbd)$Xout
for(i in 1:1000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }
> range(c(Xmat %*% w))
[1] 1.460378 3.093886        ### Narrower range

vec2=hist(c(Xmat %*% w), breaks=seq(0.85,3.45, by=0.1), prob=T)$counts
 > vec2
 [1]   0   0   0   0   0   0   5  38  48  57  38  78  75  85  87 108  88  74  75
[20]  67  43  29   5   0   0   0

### Different shape !! But it is hard to know what the right shape is:
### let's create  a matrix and then print 10 different simulated
### ouptput counts for blocks of 1000, separated by `pre-iterated'
### blocks of 5000. We print the middle 14 categories corresponding
### to values 1.45--2.85.

CountMat = array(0,c(10,26))
RangMat = array(0,c(10,2))
CountMat[1,]=vec1; CountMat[2,]=vec2; 
RangMat[1,] =c(1.293272, 3.120429)
RangMat[2,] =c(1.276512, 2.739910)

for(j in 1:8) {
  X1 = Xmat[1000,]
  for(i in 1:5000) X1= Xtrans(X1,carr,lowbd)$Xout
  for(i in 1:1000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }
  proj = c(Xmat %*% w)
  RangMat[j+2,] = range(proj)
  CountMat[j+2,] = hist(proj, breaks=seq(0.85,3.45, by=0.1), prob=T)$counts }

round(RangMat,4)
        [,1]   [,2]
 [1,] 1.2933 3.1204
 [2,] 1.2765 2.7399
 [3,] 1.6365 3.1555
 [4,] 0.8659 3.1328
 [5,] 1.5188 3.2173
 [6,] 1.2595 3.1001
 [7,] 1.4764 3.0157
 [8,] 1.5336 3.1443
 [9,] 1.2628 3.1480
[10,] 1.4747 3.2455       ### Nothing too obvious except that the
 ###  variation among the first 5 is much more than among the second 5

> c(var(RangMat[1:5,1]),var(RangMat[6:10,1]))
[1] 0.08713050 0.01695485
> c(var(RangMat[1:5,2]),var(RangMat[6:10,2]))
[1] 0.036102164 0.006957089

CountMat[,7:20]

  17  33   86   97   80   86  145  101   73  76  59  51  39  22
  5   38   48   57   38   78   75   85   87 108  88  74  75  67
  0    2   28   18   43   27   58   83  120  99  91  96  85  97
 46   38   19   38   34   51   48   69  104 126 132  68  53  27
  5    6   15   19   51   62   78  129  105  77  86  97 103  82
 21   29   29   45   61   88   96   84   54  96  89  69  94  62
 13   22   49  101   86  116   92   84  107  59  86  98  27  40
  3   15   29   21   34   41   50   79   99 118 146 125 103  65
 67   46   37   75  128  122   88   71   59  56  59  55  44  24
 16   11   16   25   37   84   81   70  103 120 142 100  67  55
     

### Does not tell so much except that there is a lot of variability!

#### >>>>> #### (Second Check)
### Now re-do the whole thing, but overlay density estimates, with a
###  different w  for  6  different  blocks of 2000  separated by 3000

w = runif(7); w = w/sqrt(sum(w*w))
> round(w,4)
[1] 0.2287 0.4977 0.4973 0.3396 0.5524 0.1187 0.1346
> range(Xmat%*%w)
[1] 1.710972 3.463750

plot(c(Xmat%*%w), type="n", xlab="Projection on w", ylab="Density",
   xlim=c(1.7,3.5), ylim=c(0,1.2), main=paste("Overlaid Density", 
   " Estimates \n 6 Blocks of 2000 sep by 3000",sep=""))

Xmat = array(1, c(2000,7))
for(j in 1:6) {
   X1 = Xmat[2000,]
   for(i in 1:3000) X1= Xtrans(X1,carr,lowbd)$Xout
   for(i in 1:2000) {
       X1 = Xtrans(X1,carr,lowbd)$Xout
       Xmat[i,] = X1 }
   proj=c(Xmat %*% w)
   tmpdns = density(proj)
   lines(tmpdns$x,tmpdns$y, lty=j) }
legend(locator(), legend=paste("Blk",1:6), lty=1:6)
### Looks OK, cannot easily see glaring differences 

#### >>>>> #### (Third Check) not done here

Overlay  acf or spectral (method="ar") plots for separate blocks of
MCMC output converted into scalar series, to see that the blocks are
distributed more or less trhe same ,just as we did above for 
densities.


#### >>>>> #### (Fourth Check, confirming uniformity on a box)

## Find a box around rep(1,7) contained in the region. If it is not
   too small, then:

### Divide this box into 2^7 = 128 pieces by conditions that each
    coord is <1 or > 1, using "sign" function

### Now for each of a very large number of MCMC generated vectors
####   (1) check that all coord's within .074 of 1, discard if not.
####   (2) for those within range, create vector of signs,then takes
###        2*s-1 is each coord and inner product with 2^(0:6) and add
###        1. This gives numerical category 1:128 which should all be
###         equally likely.

#### This calculation will do two things:
####      (i)  check uniformity within the box
####      (ii) if we compare the total number in the box compared to
####           the total number N simulated, we get an estimate of 
####           the volume of the region, whose reciprocal is the
####           normalizing constant for the density !

Counts = numeric(128)
a0 = 2^(0:6)
## assume the box half-lengths (not nec equal) are given in vector "box"

## We start with the last MCMC iterate X1:  then

for(i in 1:1e5) {
        X1= Xtrans(X1,carr,lowbd)$Xout
        if(max(abs(X1-1))<box) {
           m = sum((sign(X1-1)+1)*a0)/2 + 1
           Counts[m] = Counts[m]+1 } } 

### For additional details, see the related text-file 
###    HW18slnaux.txt  in this Directory.