Estimate normal distribution parameters from smallest N samples I have a bunch of small datasets (billions of sets of 7 samples). Each dataset represents the smallest 7 samples of a larger set of 15 values which are normally distributed. Given just the smallest 7 samples, how could I calculate the mean/stdev for the full 15 sample set efficiently?
I suspect this involves using order statistics and the beta distribution but this is going beyond my usual depth in stats.

You can do almost as well as Maximum Likelihood Estimation (MLE) using Regression on Order Statistics (RoS).  The latter is simpler to program (requiring only sorting and straightforward arithmetic) and it's thirty times faster in execution (when optimized as described below).

The idea of RoS is to fit a line--using any method you like--to the observed data on a censored probability plot.  In this image of 15 points only the lowest 7 (black) points were used to fit that line.  The plotting points are quantiles of the standard Normal distribution and (therefore) the line's intercept estimates the median (which equals the mean) and its slope estimates the standard deviation.

You can precompute the plotting points because they depend only on the common size of the datasets, $15.$  A least-squares calculation of the line's parameters is ultimately a dot product: see the code for ROS below.  You can speed it up by precomputing the variance of the plotting points and simplifying the algebra.  It is massively parallelizable.
How well does this perform?  Here are the results of applying ROS along with censored Maximum Likelihood (requiring numerical optimization) for 5,000 randomly generated Normal datasets.

The plots are similar.  Both show the estimates tend to be near the true values (plotted as the red points) and, as one would expect, positively correlated.  But the MLE results are a little bit more precise.  Here is a comparison of the standard deviations in this simulation:
    Intercept Slope
RoS     0.366 0.333
MLE     0.336 0.286

That is, RoS with 7 points is about as good as MLE with 6 points.
In any given instance, RoS and MLE tend to agree, as these scatterplots of their estimates in the simulation attest.

The bottom line is

If greater efficiency is worth the amount of information in one data point (or more), then RoS is a good choice.

Here are my R implementations of these methods and the simulation.
#
# `x` is the bottom `k` order statistics and `pp` are their corresponding
# Normal plotting points.  Returns an (intercept, slope) vector of estimates.
#
ROS <- function(x, pp) {
  s <- cov(x, pp) / var(pp)
  m <- mean(x - s * pp)
  c(Intercept = m, Slope = s)
}
#
# `x` is the bottom `k` data points and `n` is the total size of the dataset 
# they are from.   Returns an (intercept, slope) vector of estimates.
#
MLE <- function(x, n) {
  y <- sort(x, decreasing = TRUE)
  k <- length(x)
  lambda <- function(theta) { # Negative log likelihood
    mu <- theta[1]
    sigma <- exp(theta[2])
    -sum(dnorm(y, mu, sigma, log = TRUE)) - 
      (n - k) * pnorm(y[1], mu, sigma, log = TRUE, lower.tail = FALSE)
  }
  beta <- optim(c(0, 0), lambda)$par
  c(Intercept = beta[1], Slope = exp(beta[2]))
}
#
# Apply both methods to simulated datasets.
#
set.seed(17)
k <- 7
n <- 15
pp <- qnorm(seq(0, 1, length.out = 2*n+1)[2 * seq_len(n)])
sim <- replicate(5e3, {
  x <- sort(rnorm(n))        # A random sample
  y <- x[seq_len(k)]         # Its `k` lowest elements
  b <- ROS(y, pp[seq_len(k)])
  b.MLE <- MLE(y, n)
  c(b, b.MLE)                # The four estimates (in each column).
})
#
# Print the standard deviations of the four estimates.
#
(matrix(signif(apply(sim, 1, sd), 3), 2, 
 dimnames = list(c("RoS", "MLE"), c("Intercept", "Slope"))))

