Given the density function of a certain population+
$$f_\theta(x)=\frac{x+1}{\theta(\theta+1)}e^{-\frac{x}{\theta}} \qquad \text{for } x>0, \ \theta>0$$
Since $f_\theta(x)$ is part of the exponential k-parametric family, from its factorization one immediately knows that $\bar X$ is a minimal sufficient statistic. Moreover, it is complete; this is a fact you may find in several statistical inference books.
However, I cannot use this result freely in my statistics course, and I may prove by definition that $\bar X$ is complete. Here is my attempt:
If one wants to prove that if $E[f(\bar X)]=0$, then one must have $f=0$ almost everywhere, first one has to explicitly give the density of $\bar X$ to compute the expectancy of $f(\bar X)$. However, it is not trivial what this density may be, as $f_\theta(x)$ does not correspond to the density of any of the well-known continuous distributions for which $\bar X$ is easily obtainable.
So, I may define the following change of variables:
$$(Y_1,Y_2,\ldots,Y_{n-1},Y_n)=(X_1,X_2,\ldots,X_{n-1},\bar X)$$
Which has as inverse function
$$(X_1,X_2,\ldots,X_n)=(Y_1,Y_2,\ldots,nY_n-(Y_1+Y_2+\cdots+Y_{n-1}))$$
And therefore, its jacobian matrix is
$$J=\begin{bmatrix} 1 & 0 & 0 & \dots & -1 \\ 0 & 1 & 0 & \dots & -1 \\ \vdots & \vdots & \vdots & & \vdots \\ 0 & 0 & 0 & \dots & n \end{bmatrix}$$
Which gives $|\det(J)|=n$. Therefore, one has
$$f_{(Y_1,\ldots,Y_n)}(y_1,\ldots,y_n)=nf_{(X_1,\ldots,X_n)}(y_1,\ldots,y_{n-1},ny_n-(y_1+\cdots+y_{n-1}))$$
where the latter equality is nothing but the joint sample density $f_\theta(x_1,\ldots,x_n)$ at the point $(y_1,\ldots,y_{n-1},ny_n-(y_1+\cdots+y_{n-1}))$
Subtituding one has
\begin{align} & f_{(Y_1,\ldots,Y_n)}(y_1,\ldots,y_n) \\ = {} & n\frac{(y_1+1)(y_2+1)\cdots(ny_n-(y_1+y_2+\cdots+y_{n-1}))}{\theta^n(\theta+1)^n}e^\frac{-(y_1+\cdots+y_{n-1}+ny_n-y_1-\cdots-y_{n-1})}{\theta} \\ = {} & n\frac{(y_1+1)(y_2+1)\cdots(ny_n-(y_1+y_2+\cdots+y_{n-1}))}{\theta^n(\theta+1)^n} e^{-\frac{ny_n}{\theta}} \end{align}
To finally obtain $f_{\bar X}(\bar X)$ one must calculate
$$f_{\bar X}(\bar X)=\int_{C_{Y_1}}\int_{C_{Y_2}}\cdots\int_{C_{Y_{n-1}}} f_{(Y_1,\ldots,Y_n)}(y_1,\ldots,y_n) \ dy_1 \, dy_2 \ldots dy_{n-1}$$
Which can be rewritten as
$$\frac{n}{\theta^n(\theta^n+1)}e^{-n\frac{y_n}{\theta}} \int_0^\infty \int_0^\infty\cdots\int_0^\infty(y_1+1)(y_2+1)\cdots(ny_n-y_1-\cdots-y_n-1) \ dy_1 \, dy_2\cdots dy_n$$
However, I don't know how to compute that integral, nor if my previous calculations are correct.
I think this density can also be found via obtaining the characteristic function, but I am not sure if that method is any better than finding the marginal distribution as I have just tried.
After all that one should still compute $\operatorname E[f(\bar X)]$; I hope its an easier task.
Thanks in advance for your collaboration.
Note that given PDF of $X_i$ is a mixture of Exponential distribution with mean $\theta$ and Gamma($2$, $\theta$) distribution with mixture weights $\frac{1}{\theta+1}$ and $\frac{\theta}{\theta+1}$.
Introduce $n$ Bernoulli($\frac{1}{\theta+1}$) random variables $\phi_1,\ldots, \phi_n$, introduce also $n$ r.v.'s $Y_1,\ldots, Y_n$ which are Exponentially distributed with mean $\theta$, and finally introduce $n$ r.v.'s $Z_1,\ldots, Z_n$ from Gamma($2$, $\theta$) distribution. Assume all r.v.'s are mutually independent. Then $X_i$ are distributed as $$ \phi_i Y_i + (1-\phi_i)Z_i $$ and also $n\overline X=\sum_{i=1}^n X_i$ is distributed as $$ n\overline X =_d \sum_{i=1}^n (\phi_i Y_i + (1-\phi_i)Z_i) =_d S^{(1)}_\nu + S^{(2)}_{n-\nu} $$ where $\nu=\sum_{i=1}^n \phi_i$ is Binomial($n$, $\frac{1}{\theta+1}$), and two independent sums $S^{(1)}_k=Y_1+\ldots+Y_k$ and $S^{(2)}_m=Z_1+\ldots+Z_m$ come from Gamma($k,\theta$) and Gamma($2m,\theta$) distributions correspondingly.
Note that summands in the equation above are dependent through $\nu$. Note also, that given $\nu=k$, we can find conditional distribution of $n\overline X$: it is Gamma($k+2(n-k),\theta$) $=$ Gamma($2n-k,\theta$). So, the conditional PDF of $n\overline X$ given $\nu=k$ is $$ f_{n\overline X | \nu}(x|k) = \frac{x^{2n-k-1}}{\theta^{2n-k}\Gamma(2n-k)}e^{-\frac{x}{\theta}} $$ Therefore for any Borel $f(x)$ $$ \mathsf E[f(n\overline X)] = \mathsf E\left[\mathsf E\left(f(n\overline X)~\big|~\nu\right)\right]$$ $$= \sum_{k=1}^n \int_0^\infty f(x)\frac{x^{2n-k-1}}{\theta^{2n-k}\Gamma(2n-k)}e^{-\frac{x}{\theta}}\,dx\cdot \binom{n}{k}\left(\frac1{\theta+1}\right)^k \left(\frac{\theta}{\theta+1}\right)^{n-k} $$
From here you can start to prove completeness.