The following code produces the plot shown below (with a warning message saying that "byt" is not a graphical parameter). The range of forecasts are noticeably different from -200 to 300. From here, my only suggestion is that maybe there is an error elsewhere in your code.

I understand that you've got more data points, so the forecasts here will be different for the full dataset. With this in mind, your plot should still look something similar to what's shown below. If it differs significantly, check for errors elsewhere in your code.

```
# Load forecast package
library(forecast)
# Store data as a numeric vector
x <- c(568,485,360,523,514,370,332,841,685,719,647,615,389,367,731,721,819,662,581,436,394,865,805,952,759,677,429,424,598,868,888,849,707,458,426,950,806,826,804,730,541,439,1070,770,989,863,737,525,461,1041,863,989,833,783,506,496)
# Convert data to a 'ts' object
x <- ts(x, frequency = 7)
# Plot forecasts from Holt-Winters' additive method
plot(hw(x, 6), byt = "1")
```

*This answer probably should have been a comment, but I decided to make it an answer because it contains code and a plot.*

Note that the probability mass function is

$$
f(x) = P(X = x) =
\begin{cases}
\frac{16}{81}&\text{if } x=0\\
\frac{32}{81}&\text{if } x=1\\
\frac{24}{81}&\text{if } x=2\\
\frac{8}{81}&\text{if } x=3\\
\frac{1}{81}&\text{if } x=4,
\end{cases}
$$
and the cumulative distribution function is

$$
F(x) =
\begin{cases}
0 & \text{if } x<0\\
\frac{16}{81} &\text{if } 0\leq x < 1\\
\frac{48}{81} &\text{if } 1\leq x < 2\\
\frac{72}{81} &\text{if } 2\leq x < 3\\
\frac{80}{81} &\text{if } 3\leq x < 4\\
1 &\text{if } 4\leq x.
\end{cases}
$$

The quantile function $Q$ is thus

$$
Q(p) =
\begin{cases}
0 & \text{if } 0< p \leq \frac{16}{81}\\
\cdots & \cdots\\
4 & \text{if } \frac{80}{81}< p.
\end{cases}
$$

I leave the rest, e.g. the $\cdots$, to you.

**R code**

```
plot(1, xlim=c(0,1), ylim=c(-1,5), type="n",
xlab= "p", ylab= "Q(p)", main = "X ~ Bin(4,1/3)")
abline(v=c(0,1), lwd=1, lty=2, col="lightgray")
segments(x0=0,x1=16/81,y0=0,y1=0, lwd=2)
points(x=16/81,y=0, pch=20, cex=1.4)
#+++++++
segments(x0=16/81,x1=48/81,y0=1,y1=1, lwd=2)
points(x=48/81,y=1, pch=20, cex=1.4)
#+++++++
segments(x0=48/81,x1=72/81,y0=2,y1=2, lwd=2)
points(x=72/81,y=2, pch=20, cex=1.4)
#+++++++
segments(x0=72/81,x1=80/81,y0=3,y1=3, lwd=2)
points(x=80/81,y=3, pch=20, cex=1.4)
#+++++++
segments(x0=80/81,x1=81/81,y0=4,y1=4, lwd=2)
```

## Best Answer

For brevity, let $x_0$ be the adult-age and let $x$ be the age (I will assume that the age is non-negative and adult-age is positive). Your function of interest $F: \mathbb{R}_{0+} \rightarrow \mathbb{R}$ is given by:

$$\begin{align} F(x) &= \int \limits_{x_0}^x \frac{ds}{\min(s,x_0)+1} \\[6pt] &= \begin{cases} \log(x+1) - \log(x_0+1) & & \text{for } 0 \leqslant x \leqslant x_0, \\[6pt] (x-x_0)/(x_0+1) & & \text{for } x > x_0. \\[6pt] \end{cases} \end{align}$$

This function has first derivative $F'(x) = 1/(\min(x,x_0)+1)$ so it is strictly increasing and has an inverse. (It is also simple to show that this is a continuous function, so its inverse is also continuous.) For this part I will use $r=F(x)$ to represent the age-calibration. Solving for $x$ using some simple algebra on each part gives the corresponding inverse function $F^{-1}: \mathbb{R} \rightarrow \mathbb{R}_{0+}$, which is:

$$F^{-1}(r) = \begin{cases} (x_0+1) \exp(r) - 1 & & \text{for } r \leqslant 0, \\[6pt] x_0 + r(x_0+1) & & \text{for } r > 0. \\[6pt] \end{cases}$$

Computation in R:It would be preferable to program these functions without hard-coding the adult-age you wish to use. If you really want to get fancy then you can also add the derivatives to the functions as`gradient`

and`hessian`

attributes, which is useful if you use these functions later in optimisation work. You can program these functions in`R`

as follows: