Predictive distributions for data from common distributions I have data that I assume comes from some distribution, such as normal. I don't care about estimating the parameters of this distribution. Instead, I'd like to know the distribution that future observations will follow.
I think this is called a "predictive distribution". When I look it up, people write about Bayesian inference. I am not asking about Bayesian specifically, though nothing against it either. I would think this is something that frequentists use as well. What is the usual frequentist practice? Do they just pretend that the predictive distribution is the same as the original distribution with parameters set equal to the best point estimates?
I'm looking for an easy set of formulas: if the data are assumed to be from this distribution (especially normal), then this is the predictive distribution, and here is how to calculate its parameters. I'm not looking for proofs. Is there an R package for this?

Replacing parameter values with their best estimates produces a plug-in estimated sampling distribution.  This is your estimate of the true data generative process.  This can be used for producing point estimates of population percentiles.  A tolerance interval is simply a confidence interval for a population percentile.  Think of the confidence intervals that come from a time-to-event analysis for population percentiles.
In addition to estimating population percentiles you may also be interested in making a prediction about a randomly selected future observation.  This requires predictive p-values and prediction intervals.  I find tolerance intervals to be more relevant in most applications since we are often interested in performing inference on the target population, the true data generative process.  I also find it a little easier to convey percentiles and tolerance intervals to non-statisticians.  Since your question asks specifically about "prediction," below is a quote on prediction intervals from my manuscript Tolerance and Prediction Intervals for Non-normal Models.
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 constructing tolerance and prediction intervals in non-normal models.
Below is a figure from my manuscript that works well to discern between confidence intervals, tolerance intervals, and prediction intervals in a non-normal setting.  The large red density is the plug-in estimated sampling distribution, an estimate of the target population.  The dark blue confidence curve and confidence interval show p-values and confidence intervals of all levels for inference on the population mean.  The slightly wider orange prediction interval shows predictive inference on the sample mean for a future experiment.  The yellow tolerance intervals show inference on the $2.5^{th}$ and $97.5^{th}$ population percentiles.  The slightly wider orange prediction intervals show predictive inference on the sample percentiles for a future experiment.  The green $95\%$ prediction interval shows predictive inference on a single subject sampled from the target population.

