Unlike the other answerers, I think your question is perfectly clear. You are asking about the joint distribution of an iid $N(\mu,\sigma^2)$ sample $x_1,x_2,\dots,x_n$ conditional on the order statistics $x_{(1)},x_{(n)}$. Since you don't know $\mu$ and $\sigma$ it seems reasonable to use non-informative, improper location and scale priors on $\mu$ and $\sigma$, that is, $\pi(\mu,\sigma)\propto 1/\sigma$. You can certainly sample from the marginal posterior distribution of the sample (marginalising out $\mu$ and $\sigma$) under these assumptions using MCMC.
It turns out, however, that the resulting posterior distribution is identical to the distribution generated by simply simulating iid standard normal samples followed by a suitable shift and rescaling of each sample such the sample mininum and maximum matches the observed values.
The following numerical example provide evidence that the two methods indeed are equivalent. A proof follows after the example.
par(mfrow=c(1,2))
n <- 10
min = 0
max = 1
# Method 1
logposterior <- function(par, min, max, n) {
mu <- par[1]
sigma <- exp(par[2])
dnorm(min, mu, sigma, log=TRUE) +
dnorm(max, mu, sigma, log=TRUE) +
(n-2)*log(pnorm(max, mu, sigma) - pnorm(min, mu, sigma))
}
chain <- MCMCpack::MCMCmetrop1R(logposterior, theta.init=c(0,0), n=n, min=min, max=max)
#>
#>
#> @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#> The Metropolis acceptance rate was 0.57224
#> @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
posteriorsamples <- apply(chain, 1, function(par) {
mu <- par[1]
sigma <- exp(par[2])
c(0,truncnorm::rtruncnorm(n - 2, min, max, mu, sigma),1)
})
hist(posteriorsamples)
# Method 2
rescaledsamples <- replicate(1e+4, {
x <- rnorm(n)
x <- (x - min(x))/(max(x) - min(x))
x <- min + (max - min)*x
})
hist(rescaledsamples)
Created on 2024-08-15 with reprex v2.1.0
The second method can be seen as a form of fiducial inference based on the idea that the data has been generated from a standard normal iid sample $z_1,z_2,\dots,z_n$ via the transformation $x_i=\mu + \sigma z_i$. Given $x_{(1)}$ and $x_{(2)}$ and knowing the joint density of the orderstatistics $z_{(1)},z_{(2)}$, the fiducial distribution of $\mu,\sigma$ is obtained by solving
\begin{align}
x_{(1)} &= \mu + \sigma z_{(1)} \\
x_{(n)} &= \mu + \sigma z_{(n)}
\end{align}
for $\mu,\sigma$ which gives
\begin{align}
\mu &= \frac{x_{(1)}z_{(n)}-x_{(n)}z_{(1)}}{z_{(n)}-z_{(1)}} \\
\sigma &= \frac{x_{(n)}-x_{(1)}}{z_{(n)}-z_{(1)}}
\end{align}
Thinking of $z_{(1)},z_{(n)}$ as random and the observations $x_{(1)},x_{(n)}$ as fixed,
the fiducial density of $\mu,\sigma$ is then given by the usual transformation formula for joint densities giving
\begin{align}
f(\mu,\sigma)&=f(z_{(1)},z_{(n)})\left|\frac{\partial(z_{(1)},z_{(n)}))}{\partial(\mu,\sigma)}\right|
\\&\propto
\phi(z_{(1)})\phi(z_{(n)})[\Phi(z_{(n)})-\Phi(z_{(1)}]^{n-2}
\left|\begin{matrix}
-\frac1\sigma & -\frac1{\sigma^2}(x_{(1)}-\mu) \\
-\frac1\sigma & -\frac1{\sigma^2}(x_{(n)}-\mu)
\end{matrix}\right|
\\&=
\phi(z_{(1)})\phi(z_{(n)})(\Phi(z_{(n)})-\Phi(z_{(1)})^{n-2}
\frac1{\sigma^3}(x_{(n)}-x_{(1)})
\\&\propto\underbrace{\frac1{\sigma}}_{\text{Prior}\\\pi(\mu,\sigma)}
\underbrace{\frac1\sigma\phi\left(\frac{x_{(1)} - \mu}\sigma\right)\frac1\sigma\phi\left(\frac{x_{(n)} - \mu}\sigma\right)\left[\Phi(\frac{x_{(n)} - \mu}\sigma)-\Phi(\frac{x_{(1)} - \mu}\sigma)\right]^{n-2}}_{\text{Likelihood }L(\mu,\sigma)=f(x_{(1)},x_{(n)}|\mu,\sigma)}.
\end{align}
The fiducial density of $\mu,\sigma$ thus happen to be identical to the posterior density of $\mu,\sigma$ based on the above non-informative prior.