Archive for unbiased estimation

a neat (theoretical) Monte Carlo result

Posted in Books, Statistics, University life with tags , , , , on December 19, 2014 by xi'an

Mark Huber just arXived a short paper where he develops a Monte Carlo approach that bounds the probability of large errors

\mathbb{P}(|\hat\mu_t-\mu|>\epsilon\mu) < 1/\delta

by computing a lower bound on the sample size r and I wondered at the presence of μ in the bound as it indicates the approach is not translation invariant. One reason is that the standard deviation of the simulated random variables is bounded by cμ. Another reason is that Mark uses as its estimator the median

\text{med}(S_1R_1,\ldots,S_tR_t)

where the S’s are partial averages of sufficient length and the R’s are independent uniforms over (1-ε,1+ε): using those uniforms may improve the coverage of given intervals but it also means that the absolute scale of the error is multiplied by the scale of S, namely μ. I first thought that some a posteriori recentering could improve the bound but since this does not impact the variance of the simulated random variables, I doubt it is possible.

Bangalore workshop [ಬೆಂಗಳೂರು ಕಾರ್ಯಾಗಾರ]

Posted in pictures, R, Running, Statistics, Travel, University life, Wines with tags , , , , , , on July 31, 2014 by xi'an

mathdeptSecond day at the Indo-French Centre for Applied Mathematics and the workshop. Maybe not the most exciting day in terms of talks (as I missed the first two plenary sessions by (a) oversleeping and (b) running across the campus!). However I had a neat talk with another conference participant that led to [what I think are] interesting questions… (And a very good meal in a local restaurant as the guest house had not booked me for dinner!)

To wit: given a target like

\lambda \exp(-\lambda) \prod_{i=1}^n \dfrac{1-\exp(-\lambda y_i)}{\lambda}\quad (*)

the simulation of λ can be demarginalised into the simulation of

\pi (\lambda,\mathbf{z})\propto \lambda \exp(-\lambda) \prod_{i=1}^n \exp(-\lambda z_i) \mathbb{I}(z_i\le y_i)

where z is a latent (and artificial) variable. This means a Gibbs sampler simulating λ given z and z given λ can produce an outcome from the target (*). Interestingly, another completion is to consider that the zi‘s are U(0,yi) and to see the quantity

\pi(\lambda,\mathbf{z}) \propto \lambda \exp(-\lambda) \prod_{i=1}^n \exp(-\lambda z_i) \mathbb{I}(z_i\le y_i)

as an unbiased estimator of the target. What’s quite intriguing is that the quantity remains the same but with different motivations: (a) demarginalisation versus unbiasedness and (b) zi ∼ Exp(λ) versus zi ∼ U(0,yi). The stationary is the same, as shown by the graph below, the core distributions are [formally] the same, … but the reasoning deeply differs.

twoversions

Obviously, since unbiased estimators of the likelihood can be justified by auxiliary variable arguments, this is not in fine a big surprise. Still, I had not thought of the analogy between demarginalisation and unbiased likelihood estimation previously. Continue reading

trying to speed up Metropolis… and failing!

Posted in R, Statistics, University life with tags , , , on June 13, 2014 by xi'an

A while ago (but still after Iceland since I used the thorn rune as a math symbol!), I wrote the following post draft as a memo. Now that Marco Banterle, Clara Grazian and myself have completed our delayed acceptance paper, it may be of interest to some readers to see how a first attempt proved fruitless.

In the past days, I tried to speed up my student Clara’s code, reducing the number of prior evaluations in a Metropolis-Hastings algorithm by trying to reject proposals with low likelihoods before computing the corresponding prior. (This is one of those rare problems where the prior is the culprit.) My first idea was to start from the case when the new likelihood þ(θ’) was lower than the previous likelihood þ(θ)

þ(θ’) < þ(θ)

as the indicator

ℑ{u≤þ(θ’)/þ(θ)}

is an unbiased estimator of the ratio þ(θ’)/þ(θ) . When u is generated from a uniform U(0,1). (All u’s in this page will be uniform U(0,1), unless otherwise indicated.) Does this mean low values of þ(θ’) can be rejected prior to computing π(θ’)? Unfortunately no, since

ℑ{u≤þ(θ’)/þ(θ)} π(θ’)/ π(θ) ∧ 1

is not an unbiased estimator of

{þ(θ’)/þ(θ)} {π(θ’)/ π(θ)} ∧ 1

Then I considered

{þ(θ’)/þ(θ)} {π(θ’)/ π(θ)} ∧ þ(θ)/þ(θ’)}

which can be associated with the unbiased estimator

ℑ{u≤þ(θ’)/þ(θ)} {π(θ’)/ π(θ) ∧ þ(θ)/þ(θ’)}

Once more unfortunately, since

{π(θ’)/ π(θ) ∧ þ(θ)/þ(θ’)}

is not easily bounded, unless π(θ’)/π(θ) itself is bounded (a self-defeating condition!), generating an unbiased estimator of this second term is not obvious.

delayedAnother attempt was to consider an exact approximation à la Andrieu and Roberts (2009, Annals of Stat.), using an unbiased estimator of þ(θ’), namely

ℑ{þ(θ)u≤þ(θ’)} þ(θ) ℑ{þ(θ)>þ(θ’)} + þ(θ’) ℑ{þ(θ)<þ(θ’)}

and to plug this unbiased estimator in the Metropolis-Hastings acceptance ratio,

π(θ’) [ℑ{þ(θ)u’≤þ(θ’)}þ(θ)ℑ{þ(θ)>þ(θ’)} + þ(θ’)ℑ{þ(θ)<þ(θ’)}] / π(θ) [ℑ{þ(θ’)u≤þ(θ)}þ(θ’)ℑ{þ(θ’)>þ(θ)} + þ(θ)ℑ{þ(θ)>þ(θ’)}] ∧ 1

which leads to

π(θ’) [ℑ{þ(θ)u’≤þ(θ’)}þ(θ)] / π(θ) þ(θ) ∧ 1 = π(θ’) ℑ{þ(θ)u’≤þ(θ’)} / π(θ) ∧ 1

when þ(θ)>þ(θ’) and to

π(θ’) þ(θ’) / π(θ) ℑ{þ(θ’)u≤þ(θ)}þ(θ’) ∧ 1 = π(θ’) / π(θ) ℑ{þ(θ’)u≤þ(θ)} ∧ 1

otherwise. But this is not correct either. As shown by the above histogram versus target.

#data  x~B(N,p)
#posterior p|x~Be(x+7.5,N+.5-x)

T=10^4
N=100;y=32
a=7.5;b=.5

chk=NULL
p=rep(.2,T)
for (t in 2:T){

  pp=p[t-1]+runif(1,-.1,.1) #random walk
  p[t]=p[t-1] #by default rejection

  if ((pp>0)&(pp<1)){

   ratio=dbinom(y,N,pp)/dbinom(y,N,p[t-1])
   if (ratio<1){

     if ((runif(1)<ratio)&(runif(1)<
        dbeta(pp,a,b)/dbeta(p[t-1],a,b)))
                p[t]=pp

      else{

       if ((runif(1)*ratio>1)||(runif(1)<
         dbeta(pp,a,b)/dbeta(p[t-1],a,b)))
                p[t]=pp}
    }}
Follow

Get every new post delivered to your Inbox.

Join 720 other followers