Fall 2009

1. Use bisection method (either use R function uniroot( ) or write your own)
to find root(s) of the following functions:


(a) f(x) = x exp(-x) - 1  

(b) f(x) = sin(x)   (many roots) we know where the roots are, but let us see
                    how the bi-section works. 

(c) f(x) = exp( - x^2) - 0.5     (two roots)

(d) f(x) = CDF of standard normal - 0.6

(e) f(x) = x^3 - 2x + 2,  when use newton, try initial guess x0=0, and others


2. Use Newton method to find the solution of the above 5 functions.
 Try different initial values. (here you may be able to use the simple
 function newton( ). But try to do something when iteration do not converge.)

[the above are due 9/14]

3. write an R function with input: a sample of n iid observations.
The random sample is assumed to be from a Cauchy
distribution with unknown location parameter.

Output: the MLE of the location and scale paremeter.

[the above are due 9/21]

3.5. Update of 3 to including a two sample case. call this function "cauchy.reg"
     Given Y_i and x_i, so that
     Y_i is a cauchy with location parameter theta_i = (alpha+beta x_i),
     and scale parameter sigma. [sort of a regression model]
     We want the MLE of alpha, beta and sigma. [see also #4 below]


4. Version 2 of your cauchy mle function. This time assume both
   location and scale are unknown and with a documentation
   and also include in the output the log likelihood value at the MLE.
   Somehow use loglikelihood value in your interation, make sure
   each iteration increases the log likelihood value.

[the above due 9/30]


5 For the data set faithful inside R, there are two variables: 
    "eruptions" and "waiting". We shall work with the variable = 15*eruptions + waiting.
    Fit a mixture of two normals to the data. Chose your initial parameter values
    by eyeball the plotting of a histogram: 
    In R: hist( 15*eruptions + waiting ). 5 parameters total.
    Please print out the estimates after each step, to see how they converge.

    [You may want to look at R code which does EM for the mix of two 2-dim 
     normal dist. with 11 parameters and plots: 
     http://en.wikipedia.org/wiki/File:Em_old_faithful.gif
     Your task: to simplify this code to one-dim. with 5 parameters] 
     You may delete the code related to plots.

[Due 10/9]

6  Given a sample with partially observed numbers: 1, 2+, 3, 4-, 5.  

where 2+ means the actual observation is larger than 2. 
Similarly 4- means the actual observation is less than 4.

Use (nonparametric) EM algorithm to find the distribution function estimator.
Starting from the empirical distribution based on the three 
complete observations only. (i.e. probability 1/3, 1/3, 1/3 on value
1, 3, 5 ). But try experiment with other initial F(t) and see
what happens.

[If all 5 observations were known, the estimator would be
the empirical distribution based on the 5 observations.]

The correct answer for the limit prob. dist. is: locations 1, 3, 5
with prob. mass: (5 - \sqrt 5)/10,  \sqrt 5/5, (5 - \sqrt 5)/10.

Note: The limit may not be achievable in finit steps [starting from the above
itinial F], so if you are getting closer to it, it will be OK.

Or you may verify, if you start from this F( ) then after 1 EM iteration
you end up with the same thing....so this is the limit.

(6b)  add one more observation to the above sample [2.5, 4.5]; 
      i.e. we have one more observation that is only known to be
      within the interval.

[above due 10/12] 
[here: oportunity to speedup an existing R function as project]



Next Question Due Nov 6.

(7a) Suppose f(x) is a density but we can only compute it up to a constant
i.e. we can compute Cf(x) for some constant C. 

Generate a random variable X from a density g(x), 

Given X we compute the ratio   r = Cf(X)/ [ M g(X)]. 
here M is a constant we chose so that to ensure r <=1 always.
(assume this ratio is always <=1 for any X, once we pick the M)

We reject the observation X with probability (1-r). 
[i.e. accept it with probability r.]

Repeat the above procedure to get a random sample. Show that the 
random sample has distributed as f(x).


(7b) Write an R-code that uses rejection method to generate r.v. that 
follow a density given by

f(x) = C*I[-1 < x < 1]*exp( -x^4 )

where C is the standardizing constant.



Next Question Due Nov 18

(8) Variance of Median: 

First write median as the solution to an estimating equation.
Use resampling methods to estimate the variance of the median.

For definitness, let us fix a sample by

set.seed(123)
mydata2 <- rcauchy(100)

But you should pretend you do not know the sample is from a cauchy.

 To establish as a gold standard, use simulation to find the
variance of the median, this time you are allowed to use the knowledge
that the sample is from a cauchy. (so you can generate more samples)
Run enough times that took at least 5 min on your computer.























































Some questions from past.



5. Consider the function G(x1, x2) = 100( x2 - x1^2 )^2 + (1-x1)^2

Use gradian method and Newton method (either your own or use one of the
buit-in functions) to search the minimium,
starting from  ( -1.2, 1).
Compare the convergence speed, both in terms of the number of 
iterations and the total time (this is probably not very meaningful
here, but we are learning how to use system.time( ) ).


6. Given n numbers a_1, a_2, .... a_n 
  (not all the same, some negative some positive)
   Please find the p_i's that maximize   

                           \sum log p_i

   with 2 constrains:
                         \sum p_i = 1
                         \sum a_i p_i = 0

   (a) use Lagrange multiplier to reduce the problem to:
          find lambda such that

               \sum a_i/(1+ lambda a_i) = 0

   (b) then p_i = 1/n * (1/(1+lambda a_i)  is the probability we want.

[above due 10/5]


(7)  Use stochastic approx. to find the LD50 of this binomial random
   variable, using one outcome at a time.

   i.e. all we are given is   y = rbinom(1,1, prob=exp(-x) )
   but you can chose the value x.
   Our objective is to find the value x0 (dose) such that the probability
   of success is 50%.  (You are not to use the property of exp() function,
   treat rbinom() as a black box)

   We of course knew this happens when exp(-x) = 0.5, i.e.
   x= - log(0.5). But let us use stochastic approx. and only
   the binomial r.v. outcome to find the x0.

   Possible Iteration: x(n) = x(n-1) + stepsize*C*[y(n-1) - 0.5]

[above due 10/12]


resampling problems:

(9) Suppose X1, X2 ... are iid Gamma (alpha, lambda) r.v.s.
(both parameters are unknown)
We want to estimate alpha, the shape parameter by

   (\bar x)^2 / s^2

Think of many ways you can construct (approximate) confidence intervals
(90%) and try them on this estimator. Both parametric and nonparametric
confidence interval.

For definitness, let us fix a sample by using this R commend

set.seed(123)
mydata <- rgamma(100, shape=2.2)


[above due 11/11]




(13) Given n pairs of data (Y, X). Assume a linear model 

Y_i = \alpha + \beta X_i + \epsilon_i

with \epsilon_i being assumed iid N(0,1). Assume the prior distribution for
\alpha, \beta being two independent N(0, 10). 

Use R package mcmc  to generate r.v. that follow the posterior
distribution, i.e. f( \alpha , \beta | data )

Work this out with data set DNase, where n=176, Y=density, X=log(conc).
Give an estimate of the means of the posterior density.
Use at least 100000 run.

[above due 11/30]