Let $\theta \sim Gamma(1,2)$ and $X_1,...,X_n$ iid such that $X_i|\theta \sim Poisson(\theta)$. It is asked to determine the best sample size $n^*$ such that the posteriori risk

$$L(\theta, d) = (\theta-d)^2 + log(n+1)$$

is minimum. That is, we must choose such that

$$ \int_X \int_{\Theta} L(\theta, d) f(\theta, x) d\theta dx$$

is minimum for $x = (x_1,...,x_n)$. For the quadratic loss, we know that it is minimized for the decision $$\delta(X) = E(\theta|X)$$ so I must to evaluate the expectation of the posteriori for each $n$, sum with $\log(n+1)$ and check when it is minimized.It is required to calculate the expectation via Monte Carlo.

My attempt: My approach is the following: Generate 1000 values for $\theta$ and, for each $theta_i$, generate 100 values from $X \sim Poisson(\theta_i)$, $i=1,...,1000$. Then, we calculate for each $n=1,...,100$

$$\frac{1}{n}\frac{1}{1000} \sum_{i=1}^{1000} \sum_{j=1}^{n} \theta(i)x(i,j) + log(n+1)$$

The following matlab code is supposed to do this:

function [theta,X] = minRisk()

% Generate 1000 values for theta
theta = gamrnd(1,2,1,1000);
X = [];

for i = 1 : size(theta,2)
    % Simulate sample of X|theta for each theta
    X = [X; poissrnd(theta(i),1,100)];

% Each line i of this matrix contains the values of X|theta(i) generated
% multiplied by theta(i)
for i = 1 : 1000
    M(i,:) = theta(i) * X(i, :);


for n = 1:100
   E=[E; 1/(n*1000) * sum(sum(M(:,1:n))) + log(n+1)];

plot(1:100, E)


But plotting E, it have a logarithm behavior. I don't know what I am doing wrong. If one could check my code, I would be glad. R codes are welcome too.


Explanation: The goal is to minimise $$ \mathbb{E}\left[\left(\theta-\mathbb{E}\left[\theta|X_{1:n}\right]\right)^2 \right]+\log(n+1) $$ not $$ \mathbb{E}\left[\theta\bar{X}_{n}\right]+\log(n+1) $$ which corresponds to your Monte Carlo approximation $$ \frac{1}{n}\frac{1}{1000} \sum_{i=1}^{1000} \sum_{j=1}^{n} \theta(i)x(i,j) + \log(n+1) $$

Resolution: For a Gamma Ga(1,2) prior, the posterior distribution on $\theta$ is $$ \theta|x_{1:n}\sim\text{Ga}(1+n\bar{x}_n,n+2) $$ This implies that $$\mathbb{E}\left[\theta|x_{1:n}\right]=\dfrac{1+n\bar{x}_n}{2+n}$$ and that $$\text{var}\left[\theta|x_{1:n}\right]=\dfrac{1+n\bar{x}_n}{(2+n)^2}$$ Therefore $$ \mathbb{E}\left[\left(\theta-\mathbb{E}\left[\theta|X_{1:n}\right]\right)^2 \right]=\mathbb{E}\left[\text{var}\left[\theta|X_{1:n}\right]\right]=\mathbb{E}\left[\dfrac{1+n\bar{X}_n}{(2+n)^2}\right]=\dfrac{1+n\mathbb{E}\left[\bar{X}_n\right]}{(2+n)^2} $$ that can be evaluated by Monte Carlo.

Note: If, for the experiment's sake, you do not want to use the closed form expression for the variance, you do need an extra layer in your simulations, that is, simulate $\theta$, then [one] $x$ given $\theta$, then [many] $\theta'$ given $x$.

  • $\begingroup$ What if I have a mixture of gammas instead of only one gamma? In this case, I cannot use a closed form, right? So I didn't understand how should I approach. $\endgroup$ – Giiovanna Nov 17 '15 at 21:00
  • $\begingroup$ If you have a mixture of Gamma, the posterior is also a mixture of Gamma. This is closed form. From there you should be able to find a closed form for the posterior variance. If not, you can simulate from the posterior to approximate the posterior variance for each simulated value of $X_{1:n}$. Using for instance Gibbs or Metropolis algorithms when the posterior is non-standard. $\endgroup$ – Xi'an Nov 17 '15 at 21:33
  • $\begingroup$ Thanks for your help. I will try a little more and, if I have more doubts, I talk with you again! $\endgroup$ – Giiovanna Nov 17 '15 at 21:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.