# How to do estimation, when only summary statistics are available?

This is in part motivated by the following question and the discussion following it.

Suppose the iid sample is observed, $X_i\sim F(x,\theta)$. The goal is to estimate $\theta$. But original sample is not available. What we have instead are some statistics of the sample $T_1,...,T_k$. Suppose $k$ is fixed. How do we estimate $\theta$? What would be maximum likelihood estimator in this case?

• If $T_i=f(X_i)$ for a known function $f$ then you can write down the distribution of $T_i$ and the maximum likelihood estimator is derived in the usual way. But you have not precised what are the $T_i$ ? – Stéphane Laurent Sep 21 '12 at 9:44
• I am interested in the case when $T_i=f(X_1,...,X_n)$ for known $f$. This was what I meant when I said that $T_i$ are sample statistics. – mpiktas Sep 21 '12 at 10:35
• So what is the difference between $T_i$ and $T_j$ ? – Stéphane Laurent Sep 21 '12 at 11:14
• Sorry, that should have been $f_i$, not one $f$. We have several functions $f_i$, which take as an argument entire sample. – mpiktas Sep 21 '12 at 12:43
• Isn't this what maximum entropy was designed for? – probabilityislogic Oct 3 '12 at 11:53

In this case, you can consider an ABC approximation of the likelihood (and consequently of the MLE) under the following assumption/restriction:

Assumption. The original sample size $$n$$ is known.

This is not a wild assumption given that the quality, in terms of convergence, of frequentist estimators depends on the sample size, therefore one cannot obtain arbitrarily good estimators without knowing the original sample size.

The idea is to generate a sample from the posterior distribution of $$\theta$$ and, in order to produce an approximation of the MLE, you can use an importance sampling technique as in  or to consider a uniform prior on $$\theta$$ with support on a suitable set as in .

I am going to describe the method in . First of all, let me describe the ABC sampler.

ABC Sampler

Let $$f(\cdot\vert\theta)$$ be the model that generates the sample where $$\theta \in \Theta$$ is a parameter (to be estimated), $$T$$ be a statistic (a function of the sample) and $$T_0$$ be the observed statistic, in the ABC jargon this is called a summary statistic, $$\rho$$ be a metric, $$\pi(\theta)$$ a prior distribution on $$\theta$$ and $$\epsilon>0$$ a tolerance. Then, the ABC-rejection sampler can be implemented as follows.

1. Sample $$\theta^*$$ from $$\pi(\cdot)$$.
2. Generate a sample $$\bf{x}$$ of size $$n$$ from the model $$f(\cdot\vert\theta^*)$$.
3. Compute $$T^*=T({\bf x})$$.
4. If $$\rho(T^*,T_0)<\epsilon$$, accept $$\theta^*$$ as a simulation from the posterior of $$\theta$$.

This algorithm generates an approximate sample from the posterior distribution of $$\theta$$ given $$T({\bf x})=T_0$$. Therefore, the best scenario is when the statistic $$T$$ is sufficient but other statistics can be used. For a more detailed description of this see this paper.

Now, in a general framework, if one uses a uniform prior that contains the MLE in its support, then the Maximum a posteriori (MAP) coincides with Maximum Likelihood Estimator (MLE). Therefore, if you consider an appropriate uniform prior in the ABC Sampler, then you can generate an approximate sample of a posterior distribution whose MAP coincides with the MLE. The remaining step consists of estimating this mode. This problem has been discussed in CV, for instance in "Computationally efficient estimation of multivariate mode".

A toy example

Let $$(x_1,...,x_n)$$ be a sample from a $$N(\mu,1)$$ and suppose that the only information available from this sample is $$\bar{x}=\dfrac{1}{n}\sum_{j=1}^n x_j$$. Let $$\rho$$ be the Euclidean metric in $${\mathbb R}$$ and $$\epsilon=0.001$$. The following R code shows how to obtain an approximate MLE using the methods described above using a simulated sample with $$n=100$$ and $$\mu=0$$, a sample of the posterior distribution of size $$1000$$, a uniform prior for $$\mu$$ on $$(-0.3,0.3)$$, and a kernel density estimator for the estimation of the mode of the posterior sample (MAP=MLE).

# rm(list=ls())

# Simulated data
set.seed(1)
x = rnorm(100)

# Observed statistic
T0 = mean(x)

# ABC Sampler using a uniform prior

N=1000
eps = 0.001
ABCsamp = rep(0,N)
i=1

while(i < N+1){
u = runif(1,-0.3,0.3)
t.samp = rnorm(100,u,1)
Ts = mean(t.samp)
if(abs(Ts-T0)<eps){
ABCsamp[i]=u
i=i+1
print(i)
}
}

# Approximation of the MLE
kd = density(ABCsamp)
kd$$x[which(kd$$y==max(kd$y))]  As you can see, using a small tolerance we get a very good approximation of the MLE (which in this trivial example can be calculated from the statistic given that it is sufficient). It is important to notice that the choice of the summary statistic is crucial. Quantiles are typically a good choice for the summary statistic, but not all the choices produce a good approximation. It may be the case that the summary statistic is not very informative and then the quality of the approximation might be poor, which is well-known in the ABC community. Update: A similar approach was recently published in Fan et al. (2012). See this entry for a discussion on the paper. • (+1) For stating the correct result about the relationship between MLE and MAP and for the warning in the last paragraph (among other reasons). To make that warning more explicit, this (or any!) approach will fail miserably if the statistics at hand are ancillary or nearly so. One can consider your toy example and$T = \sum_i (X_i - \bar X)^2$, for example. – cardinal Sep 21 '12 at 10:40 • +1 @procrastinator I was going to simple say yes you can use the sufficient statistics if they are available for your model. But your extensive answers seems to have covered that. – Michael R. Chernick Sep 21 '12 at 12:24 • One simple question, you mention that uniform prior must contain MLE in its support. But MLE is a random variable which is only stochastically bounded, i.e. it can be outside of any bounded set with positive probability. – mpiktas Sep 21 '12 at 12:59 • @mpiktas For a specific sample, you have to choose the appropriate support of the uniform prior. This may change if you change the sample. It is important to note that this is not a Bayesian procedure, we are just using it as a numerical method, therefore there is no problem on playing with the choice of the prior. The smaller the support of the prior, the better. This would increase the speed of the ABC sampler but when your information is vague in the sense that you do not have a reliable clue on where the MLE is located, then you might need a larger support (and will pay the price). – user10525 Sep 21 '12 at 13:06 • @mpiktas In the toy example, you can use, for instance, a uniform prior with support on$(-1000000,1000000)$or a uniform prior with support on$(0.1,0.15)$obtaining the same results but with extremely different acceptance rates. The choice of this support is ad hoc and it is impossible to come up with a general-purpose prior given that the MLE is not stochastically bounded, as you mention. This choice can be considered as a lever of the method that has to be adjusted in each particular case. – user10525 Sep 21 '12 at 13:10 It all depends on whether or not the joint distribution of those$T_i$'s is known. If it is, e.g., $$(T_1,\ldots,T_k)\sim g(t_1,\ldots,t_k|\theta,n)$$ then you can conduct maximum likelihood estimation based on this joint distribution. Note that, unless$(T_1,\ldots,T_k)$is sufficient, this will almost always be a different maximum likelihood than when using the raw data$(X_1,\ldots,X_n)$. It will necessarily be less efficient, with a larger asymptotic variance. If the above joint distribution with density$g$is not available, the solution proposed by Procrastinator is quite appropriate. The (frequentist) maximum likelihood estimator is as follows: For$F$in the exponential family, and if your statistics are sufficient your likelihood to be maximised can always be written in the form: $$l(\theta| T) = \exp\left( -\psi(\theta) + \langle T,\phi(\theta) \rangle \right),$$ where$\langle \cdot, \cdot\rangle$is the scalar product,$T$is the vector of suff. stats. and$\psi(\cdot)$and$\phi(\cdot)\$ are continuous twice-differentiable.

The way you actually maximize the likelihood depends mostly on the possiblity to write the likelihood analytically in a tractable way. If this is possible you will be able to consider general optimisation algorithms (newton-raphson, simplex...). If you do not have a tractable likelihood, you may find it easier to compute a conditional expection as in the EM algorithm, which will also yield maximum likelihood estimates under rather affordable hypotheses.

Best

• For problems I am interested in, analytical tractability is not possible. – mpiktas Oct 8 '12 at 7:03
• The reason for non-tractability then conditions the optimization scheme. However, extensions of the EM usually allow to get arround most of these reasons. I don"t think I can be more specific in my suggestions without seeing the model itself – julien stirnemann Oct 8 '12 at 10:18