You are interested in predictive p-values and prediction intervals.  Below is a quote from my manuscript *[Tolerance and Prediction Intervals for Non-normal Models.][1]*  

In repeated sampling a prediction interval covers a future observation of a random process $100(1-\alpha)\%$ of the time. For normally distributed data $\boldsymbol{Y}_n=Y_1,...,Y_n$ when the population variance $\sigma^2$ is known the pivotal quantity $(\bar{Y}_n-Y_{n+1})/\sigma\sqrt{1/n+1}$ is ancillary since it and its sampling distribution, $N(0,1)$, do not depend on the unknown mean $\mu$, where $\bar{Y}_n=(1/n)\sum_{i=1}^n Y_i$. When pivoted this quantity results in the interval estimate $\bar{y}_n\pm z_{1-\alpha/2}\cdot\sigma\sqrt{1/n+1}$, where $z_{1-\alpha/2}$ is the $100(1-\alpha/2)^{th}$ percentile of the standard normal distribution. This is a prediction interval for the as of yet unobserved $y_{n+1}$. When $\sigma^2$ is not known the ancillary pivotal quantity of choice becomes $(\bar{Y}_n-Y_{n+1})/S\sqrt{1/n+1}\sim T_{n-1}$, where $S^2$ is the bias corrected sample variance.  In repeated sampling $\bar{y}_n\pm t_{n-1,1-\alpha/2}\cdot s\sqrt{1/n+1}$ will cover the $n+1^{th}$ observation $100(1-\alpha)\%$ of the time, regardless of the unknown fixed true $\mu$ and $\sigma^2$. The p-value testing the hypothesis $H_0$: $y_{n+1} \le c$ is given by $P\big(T_{n-1}\ge (\bar{y}_n-c)/s\sqrt{1/n+1}\big)$, the probability of the difference (discrepancy) between the observed result and the hypothesized future result or something more extreme, if these share unknown fixed parameters $\mu$ and $\sigma$.  This probability forms the level of confidence that $y_{n+1}$ will be less than or equal to $c$, and is useful for controlling the type I error rate $\alpha$ when predicting $y_{n+1}$. The upper p-value function of all upper-tailed predictive p-values as a function of the hypothesis being tested is $H(y_{n+1})=1-\Phi_{n-1}\big(\hspace{1mm}\big[\bar{y}_n- y_{n+1}\big]\big/s\sqrt{1/n+1}\hspace{1mm} \big)$, where $\Phi_{n-1}$ denotes the cdf of a $T_{n-1}$ random variable. The null value $c$ is replaced with $y_{n+1}$ to denote that this is a function of all possible hypotheses around $y_{n+1}$. One can analogously define $H^-(y_{n+1})$ as the function of all lower-tailed p-values. The corresponding prediction confidence curve defined as 
\begin{eqnarray}
C(y_{n+1})&\equiv& \left\{ \begin{array}{cc}
H(y_{n+1}) & \text{if } y_{n+1}\le\bar{y}_n \\
 &  \nonumber\\
 H^{-}(y_{n+1})  & \text{if }y_{n+1}\ge\bar{y}_n, \end{array}  \right.\nonumber
\end{eqnarray}
and prediction confidence density $h(y_{n+1})\equiv dH(y_{n+1})/dy_{n+1}$, depict p-values and prediction intervals of all levels for hypotheses around $y_{n+1}$.

The above approach is easily modified if you are interested in predicting, say, $\bar{y}_{N-n}=\frac{1}{N-n}\sum_{i=n+1}^N y_i$ instead of ${y}_{n+1}$.  My paper also discusses practical approximate ancillary pivotal quantities and other methods for prediction in non-normal models.  


  [1]: https://www.researchgate.net/publication/350055270_Tolerance_and_Prediction_Intervals_for_Non-normal_Models