We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
DeepAR with NegativeBinomial cannot generate values above 1e6
## Description
A DeepAR model with NegativeBinomial output distribution cannot generate values significantly above 1e6.

## To Reproduce
I attach a jupyter notebook where I generate artificial timeseries with values between 0 and 1e8, train a model and plot the forecast. I compressed the notebook with zip as .ipynb files are not supported as attachments.

[1e6.ipynb.zip](https://github.com/awslabs/gluon-ts/files/8069187/1e6.ipynb.zip)

## Error message or code output
Please see the attached notebook.

![image](https://user-images.githubusercontent.com/16802400/154065593-3cfa3573-f5d5-4baf-a87e-714400f6e804.png)

## Environment
- Operating system: Ubuntu 20.04, linux kernel 5.13.0-28-generic
- Python version: 3.8.10
- GluonTS version: 0.8.1
- MXNet version: 1.9.0

I vaguely recall that 
I observed this issue alredy in gluonts versions 0.4.x.

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `src/gluonts/mx/distribution/neg_binomial.py`
Content:
```
1 # Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 #
3 # Licensed under the Apache License, Version 2.0 (the "License").
4 # You may not use this file except in compliance with the License.
5 # A copy of the License is located at
6 #
7 #     http://www.apache.org/licenses/LICENSE-2.0
8 #
9 # or in the "license" file accompanying this file. This file is distributed
10 # on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
11 # express or implied. See the License for the specific language governing
12 # permissions and limitations under the License.
13 
14 from typing import Dict, List, Optional, Tuple
15 
16 import numpy as np
17 
18 from gluonts.core.component import validated
19 from gluonts.mx import Tensor
20 
21 from .deterministic import DeterministicOutput
22 from .distribution import Distribution, _sample_multiple, getF, softplus
23 from .distribution_output import DistributionOutput
24 from .mixture import MixtureDistributionOutput
25 
26 
27 class NegativeBinomial(Distribution):
28     r"""
29     Negative binomial distribution, i.e. the distribution of the number of
30     successes in a sequence of independent Bernoulli trials.
31 
32     Parameters
33     ----------
34     mu
35         Tensor containing the means, of shape `(*batch_shape, *event_shape)`.
36     alpha
37         Tensor of the shape parameters, of shape `(*batch_shape, *event_shape)`.
38     F
39     """
40 
41     is_reparameterizable = False
42 
43     @validated()
44     def __init__(self, mu: Tensor, alpha: Tensor) -> None:
45         self.mu = mu
46         self.alpha = alpha
47 
48     @property
49     def F(self):
50         return getF(self.mu)
51 
52     @property
53     def batch_shape(self) -> Tuple:
54         return self.mu.shape
55 
56     @property
57     def event_shape(self) -> Tuple:
58         return ()
59 
60     @property
61     def event_dim(self) -> int:
62         return 0
63 
64     def log_prob(self, x: Tensor) -> Tensor:
65         alphaInv = 1.0 / self.alpha
66         alpha_times_mu = self.alpha * self.mu
67         F = self.F
68         ll = (
69             x * F.log(alpha_times_mu / (1.0 + alpha_times_mu))
70             - alphaInv * F.log1p(alpha_times_mu)
71             + F.gammaln(x + alphaInv)
72             - F.gammaln(x + 1.0)
73             - F.gammaln(alphaInv)
74         )
75         return ll
76 
77     @property
78     def mean(self) -> Tensor:
79         return self.mu
80 
81     @property
82     def stddev(self) -> Tensor:
83         return self.F.sqrt(self.mu * (1.0 + self.mu * self.alpha))
84 
85     def sample(
86         self, num_samples: Optional[int] = None, dtype=np.float32
87     ) -> Tensor:
88         def s(mu: Tensor, alpha: Tensor) -> Tensor:
89             F = self.F
90             tol = 1e-5
91             r = 1.0 / alpha
92             theta = alpha * mu
93             r = F.minimum(F.maximum(tol, r), 1e10)
94             theta = F.minimum(F.maximum(tol, theta), 1e10)
95             x = F.minimum(F.random.gamma(r, theta), 1e6)
96             return F.random.poisson(lam=x, dtype=dtype)
97 
98         return _sample_multiple(
99             s, mu=self.mu, alpha=self.alpha, num_samples=num_samples
100         )
101 
102     @property
103     def args(self) -> List:
104         return [self.mu, self.alpha]
105 
106 
107 class NegativeBinomialOutput(DistributionOutput):
108     args_dim: Dict[str, int] = {"mu": 1, "alpha": 1}
109     distr_cls: type = NegativeBinomial
110 
111     @classmethod
112     def domain_map(cls, F, mu, alpha):
113         epsilon = np.finfo(cls._dtype).eps  # machine epsilon
114 
115         mu = softplus(F, mu) + epsilon
116         alpha = softplus(F, alpha) + epsilon
117         return mu.squeeze(axis=-1), alpha.squeeze(axis=-1)
118 
119     # Overwrites the parent class method.
120     # We cannot scale using the affine transformation since negative binomial should return integers.
121     # Instead we scale the parameters.
122     def distribution(
123         self,
124         distr_args,
125         loc: Optional[Tensor] = None,
126         scale: Optional[Tensor] = None,
127     ) -> NegativeBinomial:
128         mu, alpha = distr_args
129         if scale is None:
130             return NegativeBinomial(mu, alpha)
131         else:
132             F = getF(mu)
133             mu = F.broadcast_mul(mu, scale)
134             return NegativeBinomial(mu, alpha, F)
135 
136     @property
137     def event_shape(self) -> Tuple:
138         return ()
139 
140 
141 def ZeroInflatedNegativeBinomialOutput() -> MixtureDistributionOutput:
142     return MixtureDistributionOutput(
143         distr_outputs=[NegativeBinomialOutput(), DeterministicOutput(0)]
144     )
145
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/src/gluonts/mx/distribution/neg_binomial.py b/src/gluonts/mx/distribution/neg_binomial.py
--- a/src/gluonts/mx/distribution/neg_binomial.py
+++ b/src/gluonts/mx/distribution/neg_binomial.py
@@ -87,13 +87,9 @@
     ) -> Tensor:
         def s(mu: Tensor, alpha: Tensor) -> Tensor:
             F = self.F
-            tol = 1e-5
             r = 1.0 / alpha
             theta = alpha * mu
-            r = F.minimum(F.maximum(tol, r), 1e10)
-            theta = F.minimum(F.maximum(tol, theta), 1e10)
-            x = F.minimum(F.random.gamma(r, theta), 1e6)
-            return F.random.poisson(lam=x, dtype=dtype)
+            return F.random.poisson(lam=F.random.gamma(r, theta), dtype=dtype)
 
         return _sample_multiple(
             s, mu=self.mu, alpha=self.alpha, num_samples=num_samples

{"golden_diff": "diff --git a/src/gluonts/mx/distribution/neg_binomial.py b/src/gluonts/mx/distribution/neg_binomial.py\n--- a/src/gluonts/mx/distribution/neg_binomial.py\n+++ b/src/gluonts/mx/distribution/neg_binomial.py\n@@ -87,13 +87,9 @@\n     ) -> Tensor:\n         def s(mu: Tensor, alpha: Tensor) -> Tensor:\n             F = self.F\n-            tol = 1e-5\n             r = 1.0 / alpha\n             theta = alpha * mu\n-            r = F.minimum(F.maximum(tol, r), 1e10)\n-            theta = F.minimum(F.maximum(tol, theta), 1e10)\n-            x = F.minimum(F.random.gamma(r, theta), 1e6)\n-            return F.random.poisson(lam=x, dtype=dtype)\n+            return F.random.poisson(lam=F.random.gamma(r, theta), dtype=dtype)\n \n         return _sample_multiple(\n             s, mu=self.mu, alpha=self.alpha, num_samples=num_samples\n", "issue": "DeepAR with NegativeBinomial cannot generate values above 1e6\n## Description\r\nA DeepAR model with NegativeBinomial output distribution cannot generate values significantly above 1e6.\r\n\r\n## To Reproduce\r\nI attach a jupyter notebook where I generate artificial timeseries with values between 0 and 1e8, train a model and plot the forecast. I compressed the notebook with zip as .ipynb files are not supported as attachments.\r\n\r\n[1e6.ipynb.zip](https://github.com/awslabs/gluon-ts/files/8069187/1e6.ipynb.zip)\r\n\r\n## Error message or code output\r\nPlease see the attached notebook.\r\n\r\n![image](https://user-images.githubusercontent.com/16802400/154065593-3cfa3573-f5d5-4baf-a87e-714400f6e804.png)\r\n\r\n## Environment\r\n- Operating system: Ubuntu 20.04, linux kernel 5.13.0-28-generic\r\n- Python version: 3.8.10\r\n- GluonTS version: 0.8.1\r\n- MXNet version: 1.9.0\r\n\r\nI vaguely recall that \r\nI observed this issue alredy in gluonts versions 0.4.x.\n", "before_files": [{"content": "# Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved.\n#\n# Licensed under the Apache License, Version 2.0 (the \"License\").\n# You may not use this file except in compliance with the License.\n# A copy of the License is located at\n#\n#     http://www.apache.org/licenses/LICENSE-2.0\n#\n# or in the \"license\" file accompanying this file. This file is distributed\n# on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either\n# express or implied. See the License for the specific language governing\n# permissions and limitations under the License.\n\nfrom typing import Dict, List, Optional, Tuple\n\nimport numpy as np\n\nfrom gluonts.core.component import validated\nfrom gluonts.mx import Tensor\n\nfrom .deterministic import DeterministicOutput\nfrom .distribution import Distribution, _sample_multiple, getF, softplus\nfrom .distribution_output import DistributionOutput\nfrom .mixture import MixtureDistributionOutput\n\n\nclass NegativeBinomial(Distribution):\n    r\"\"\"\n    Negative binomial distribution, i.e. the distribution of the number of\n    successes in a sequence of independent Bernoulli trials.\n\n    Parameters\n    ----------\n    mu\n        Tensor containing the means, of shape `(*batch_shape, *event_shape)`.\n    alpha\n        Tensor of the shape parameters, of shape `(*batch_shape, *event_shape)`.\n    F\n    \"\"\"\n\n    is_reparameterizable = False\n\n    @validated()\n    def __init__(self, mu: Tensor, alpha: Tensor) -> None:\n        self.mu = mu\n        self.alpha = alpha\n\n    @property\n    def F(self):\n        return getF(self.mu)\n\n    @property\n    def batch_shape(self) -> Tuple:\n        return self.mu.shape\n\n    @property\n    def event_shape(self) -> Tuple:\n        return ()\n\n    @property\n    def event_dim(self) -> int:\n        return 0\n\n    def log_prob(self, x: Tensor) -> Tensor:\n        alphaInv = 1.0 / self.alpha\n        alpha_times_mu = self.alpha * self.mu\n        F = self.F\n        ll = (\n            x * F.log(alpha_times_mu / (1.0 + alpha_times_mu))\n            - alphaInv * F.log1p(alpha_times_mu)\n            + F.gammaln(x + alphaInv)\n            - F.gammaln(x + 1.0)\n            - F.gammaln(alphaInv)\n        )\n        return ll\n\n    @property\n    def mean(self) -> Tensor:\n        return self.mu\n\n    @property\n    def stddev(self) -> Tensor:\n        return self.F.sqrt(self.mu * (1.0 + self.mu * self.alpha))\n\n    def sample(\n        self, num_samples: Optional[int] = None, dtype=np.float32\n    ) -> Tensor:\n        def s(mu: Tensor, alpha: Tensor) -> Tensor:\n            F = self.F\n            tol = 1e-5\n            r = 1.0 / alpha\n            theta = alpha * mu\n            r = F.minimum(F.maximum(tol, r), 1e10)\n            theta = F.minimum(F.maximum(tol, theta), 1e10)\n            x = F.minimum(F.random.gamma(r, theta), 1e6)\n            return F.random.poisson(lam=x, dtype=dtype)\n\n        return _sample_multiple(\n            s, mu=self.mu, alpha=self.alpha, num_samples=num_samples\n        )\n\n    @property\n    def args(self) -> List:\n        return [self.mu, self.alpha]\n\n\nclass NegativeBinomialOutput(DistributionOutput):\n    args_dim: Dict[str, int] = {\"mu\": 1, \"alpha\": 1}\n    distr_cls: type = NegativeBinomial\n\n    @classmethod\n    def domain_map(cls, F, mu, alpha):\n        epsilon = np.finfo(cls._dtype).eps  # machine epsilon\n\n        mu = softplus(F, mu) + epsilon\n        alpha = softplus(F, alpha) + epsilon\n        return mu.squeeze(axis=-1), alpha.squeeze(axis=-1)\n\n    # Overwrites the parent class method.\n    # We cannot scale using the affine transformation since negative binomial should return integers.\n    # Instead we scale the parameters.\n    def distribution(\n        self,\n        distr_args,\n        loc: Optional[Tensor] = None,\n        scale: Optional[Tensor] = None,\n    ) -> NegativeBinomial:\n        mu, alpha = distr_args\n        if scale is None:\n            return NegativeBinomial(mu, alpha)\n        else:\n            F = getF(mu)\n            mu = F.broadcast_mul(mu, scale)\n            return NegativeBinomial(mu, alpha, F)\n\n    @property\n    def event_shape(self) -> Tuple:\n        return ()\n\n\ndef ZeroInflatedNegativeBinomialOutput() -> MixtureDistributionOutput:\n    return MixtureDistributionOutput(\n        distr_outputs=[NegativeBinomialOutput(), DeterministicOutput(0)]\n    )\n", "path": "src/gluonts/mx/distribution/neg_binomial.py"}], "after_files": [{"content": "# Copyright 2018 Amazon.com, Inc. or its affiliates. All Rights Reserved.\n#\n# Licensed under the Apache License, Version 2.0 (the \"License\").\n# You may not use this file except in compliance with the License.\n# A copy of the License is located at\n#\n#     http://www.apache.org/licenses/LICENSE-2.0\n#\n# or in the \"license\" file accompanying this file. This file is distributed\n# on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either\n# express or implied. See the License for the specific language governing\n# permissions and limitations under the License.\n\nfrom typing import Dict, List, Optional, Tuple\n\nimport numpy as np\n\nfrom gluonts.core.component import validated\nfrom gluonts.mx import Tensor\n\nfrom .deterministic import DeterministicOutput\nfrom .distribution import Distribution, _sample_multiple, getF, softplus\nfrom .distribution_output import DistributionOutput\nfrom .mixture import MixtureDistributionOutput\n\n\nclass NegativeBinomial(Distribution):\n    r\"\"\"\n    Negative binomial distribution, i.e. the distribution of the number of\n    successes in a sequence of independent Bernoulli trials.\n\n    Parameters\n    ----------\n    mu\n        Tensor containing the means, of shape `(*batch_shape, *event_shape)`.\n    alpha\n        Tensor of the shape parameters, of shape `(*batch_shape, *event_shape)`.\n    F\n    \"\"\"\n\n    is_reparameterizable = False\n\n    @validated()\n    def __init__(self, mu: Tensor, alpha: Tensor) -> None:\n        self.mu = mu\n        self.alpha = alpha\n\n    @property\n    def F(self):\n        return getF(self.mu)\n\n    @property\n    def batch_shape(self) -> Tuple:\n        return self.mu.shape\n\n    @property\n    def event_shape(self) -> Tuple:\n        return ()\n\n    @property\n    def event_dim(self) -> int:\n        return 0\n\n    def log_prob(self, x: Tensor) -> Tensor:\n        alphaInv = 1.0 / self.alpha\n        alpha_times_mu = self.alpha * self.mu\n        F = self.F\n        ll = (\n            x * F.log(alpha_times_mu / (1.0 + alpha_times_mu))\n            - alphaInv * F.log1p(alpha_times_mu)\n            + F.gammaln(x + alphaInv)\n            - F.gammaln(x + 1.0)\n            - F.gammaln(alphaInv)\n        )\n        return ll\n\n    @property\n    def mean(self) -> Tensor:\n        return self.mu\n\n    @property\n    def stddev(self) -> Tensor:\n        return self.F.sqrt(self.mu * (1.0 + self.mu * self.alpha))\n\n    def sample(\n        self, num_samples: Optional[int] = None, dtype=np.float32\n    ) -> Tensor:\n        def s(mu: Tensor, alpha: Tensor) -> Tensor:\n            F = self.F\n            r = 1.0 / alpha\n            theta = alpha * mu\n            return F.random.poisson(lam=F.random.gamma(r, theta), dtype=dtype)\n\n        return _sample_multiple(\n            s, mu=self.mu, alpha=self.alpha, num_samples=num_samples\n        )\n\n    @property\n    def args(self) -> List:\n        return [self.mu, self.alpha]\n\n\nclass NegativeBinomialOutput(DistributionOutput):\n    args_dim: Dict[str, int] = {\"mu\": 1, \"alpha\": 1}\n    distr_cls: type = NegativeBinomial\n\n    @classmethod\n    def domain_map(cls, F, mu, alpha):\n        epsilon = np.finfo(cls._dtype).eps  # machine epsilon\n\n        mu = softplus(F, mu) + epsilon\n        alpha = softplus(F, alpha) + epsilon\n        return mu.squeeze(axis=-1), alpha.squeeze(axis=-1)\n\n    # Overwrites the parent class method.\n    # We cannot scale using the affine transformation since negative binomial should return integers.\n    # Instead we scale the parameters.\n    def distribution(\n        self,\n        distr_args,\n        loc: Optional[Tensor] = None,\n        scale: Optional[Tensor] = None,\n    ) -> NegativeBinomial:\n        mu, alpha = distr_args\n        if scale is None:\n            return NegativeBinomial(mu, alpha)\n        else:\n            F = getF(mu)\n            mu = F.broadcast_mul(mu, scale)\n            return NegativeBinomial(mu, alpha, F)\n\n    @property\n    def event_shape(self) -> Tuple:\n        return ()\n\n\ndef ZeroInflatedNegativeBinomialOutput() -> MixtureDistributionOutput:\n    return MixtureDistributionOutput(\n        distr_outputs=[NegativeBinomialOutput(), DeterministicOutput(0)]\n    )\n", "path": "src/gluonts/mx/distribution/neg_binomial.py"}]}