Determine the best sample size for minimum expected loss

Question

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.

Answer 1

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$.