Small Log Concerning Inverse Functions in R           
===========================================      2/20/08

Suppose that we take an existing function F and 
want to code its inverse as an R function. 

> F = plogis   #### we take as example the logistic df
  ### which has the "logit" log(p/(1-p)) as inverse function
> xt = runif(100,.1,.99)
  sum(abs( -log(1/xt-1) - qlogis(xt)))
[1] 2.030808e-14    
     ### This shows we have correct form of inverse

### Now let us pretend that F does NOT have an analytically
###  defined inverse, and code using uniroot:

> Finv = function(u) uniroot(function(x,.u) F(x) - .u, c(-10,10),
       tol=1e-9, .u = u)$root for(i in 1:100)
       cat(round(Finv(xt[i])-qlogis(xt[i]),8),".  ")

  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .
  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .
  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .
  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .
  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  
  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 . 
  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  0 .  
  0 .  0 .


---------------------------------------------------------------------

Additional Material on Simulating Mixtures
------------------------------------------

A. Simulating from Piecewise Uniform density

This can be viewed as a simulation from a mixture of uniform densities
with disjoint supports.

> PieceUnif = function(Nrep, avec, pvec) {
### Simulate Nrep variates from piecewise uniform density
### a_1,a_2,...,a_{n+1}  are the entries of  avec  and represent
###      the successive interval endpoints
### p_1,.., p_n  are probabilities summing to 1, for the random 
###      variable of interest to fall in the successive intervals
###      (a_i,a_{i+1}]
   inds = sample(1:length(pvec), Nrep, rep=T, prob=pvec)
   uvec = runif(Nrep)
   avec[inds] + uvec* (avec[inds+1]-avec[inds]) }

> hist(PieceUnif(1000, seq(0,1,by=.1), pvec=rep(c(.06,.14),5)), 
        nclass=20)
#### Shows the desired sawtoothed behavior !


B. Simulating from Piecwise Triangular Density
[not completely tested yet ...]

> PieceTriang = function(Nrep, avec, pvec, slop) {
### Simulate Nrep variates from piecewise linear density
### avec and pvec have the same descriptions as in PieceUnif
### slop is a vector of length n (same length as pvec)
### containing the slopes slop_i  of the linear density 
### segment over the interval (a_i,a_{i+1}].
### NOTE: piecewise linear density segment on the interval
###  (a_i,a_{i+1}] is taken proportional to
###    1/(a_{i+1}-a_i)-0.5*slop_i + 
###             slop_i*(x-a_i)/(a_{i+1}-a_i)
   inds = sample(1:length(pvec), Nrep, rep=T, prob=pvec)
   uvec = runif(Nrep)
   wid = avec[inds+1]-avec[inds]
   svec = slop[inds]
   Zvec = avec[inds]
### output is vector of Nrep simulated variates
   zind = (svec==0)
   nzind = !zind
   Zvec[zind] = Zvec[zind] + wid[zind]*uvec[zind]
   aux = sqrt(2*wid[nzind]*uvec[nzind]/svec[nzind] +
         (1/svec[nzind]-wid[nzind]/2)^2)
   Zvec[nzind] = Zvec[nzind] + wid[nzind]/2-
         1/svec[nzind] + (2*(svec[nzind] > 0)-1)*aux
   Zvec  }

TWO TESTS:

> hist(PieceTriang(1000, c(0,.5,1), c(.4,.6), c(2, 1)),
         nclass=8, xlab="Z axis",  main=
  "Simulated Variates, 2-Piece Trapedzoidal Density")
  ### histogram has clear upward slope on [0,.5] and 
  ###    gentler upward slope on (.5,1].

> hist(PieceTriang(10000, c(0,.3,.7,1), c(.2,.5,.3), 
       c(-1,3, 0)), nclass=20, xlab="Z axis",
      main="3 Piece Density -- 10000 variates")
 ## Here the leftmost piece of the density has  
##    negative slope, and the three distinct pieces 
##    can be seen clearly in the histogram.