# Gamma Regression as the Last Layer of the Neural Network

My current task involves predicting data that follows a Gamma distribution. To avoid confusion of notations, in the following discussion, the p.d.f will be $$\mathbb{P}(y|\alpha, \beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}y^{\alpha - 1}\exp(-\beta y)$$ and by taking $$\eta = -\beta, b(y) = \frac{y^{\alpha - 1}}{\Gamma(\alpha)}, a(\eta) = -\alpha \log(-\eta)$$, the activation will be $$h_\theta(x) = -\frac{\alpha}{\theta^T x}$$ and the negative log-likelihood loss will be $$l(\theta) = \alpha \log h_\theta (x) + \frac{\alpha y}{h_\theta(x)}+\log \frac{\Gamma(\alpha)}{\alpha^\alpha y^{\alpha - 1}}$$

Currently, I am using trivial activation and the traditional squared error loss in the last layer, though knowing that the distribution is by no means normal. And the model is indeed not performing optimally at the samples with the lower value, roughly around the maximum of the p.d.f. Therefore I assume that using gamma regression will help reduce the error.

If I apply ReLu activation in the last but one layer, therefore the input of the last layer $$x$$ is guaranteed to be positive (for all entries), therefore to prevent "explosion" caused by the $$\theta^T x\geq 0$$, $$\theta$$ has better be negative at all times. How do I enforce these constraints (for initialization and updates) in practice (by implementing the activation and loss function in TensorFlow, for example)? Or is there any package that provides such layers? (tfp seems to be supporting the GLM models but is there a way to combine them into a deep learning model?)