# Determine the best sample size for minimum expected loss

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)];
end

% 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, :);
end

E=[];

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

plot(1:100, E)


end

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$.
• 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. – Xi'an Nov 17 '15 at 21:33