A friend helpfully suggested that this is similar to a Brownian Bridge, and I was able to use that to derive an analytical solution.
Conditions: $w$ = min boundary, $x$ = starting value, $y$ = final value, $z$ = max boundary
1. Incorporating the x and y conditions using a Brownian Bridge
The generalized Brownian Bridge is a specific Wiener Process, defined as a time-dependent normal distribution:
$$B(t)=N(\mu_t,\sigma^2_t)$$
that is tethered at either end by the conditions:
$$B(t_1)=x,\; B(t_2)=y$$
The mean and variance for this distribution are:
$$\mu_t=x+\frac{t-t_1}{t_2-t_1}(y-x)$$
$$\sigma^2_t=\frac{(t_2-t)(t-t_1)}{t_2-t_1}$$
When our time interval is $[0,1]$, these reduce to:
$$\mu_t=x(1-t) + yt$$
$$\sigma^2_t=t(1-t)$$
But then we have to condition our distribution on the observed min and max values, $w$ and $z$.
2. Min/Max Conditional
The total distribution of points observed across the entire $t\in[0,1]$, $B_T$ will also be normal (Brownian Bridges are proven to be normally distributed at all scales). We want to find the normal distribution that would be expected to have $w$ and $z$ as its min and max observed values.
It can be proven using the joint cumulative distribution functions (CDFs) that $B_T$ has a mean $\mu_T=\frac{1}2(w+z)$ halfway between the endpoints, and a standard deviation $\sigma_T=k(n)*(z-w)$ that depends on the sample size $n$.
For example, a well-known result derived from the normal CDFs is that 99.73% of the normal probability distribution lies within $3\sigma$ on either side of the mean. If we assume that $w$ and $z$ are the only two values observed outside this range, then $0.0027 = \frac{2}n$, giving $n=741$. Since this is the case where $2*3\sigma_T=z-w$, we can use this in our previous standard deviation equation to find $k(741)=\frac{1}6$.
Generalizing, where $p$ is the tail of the normal distribution containing a single extremum, and $z_p$ is the quantile:
$$2(1-p) = \frac{2}n$$
$$p=\frac{n-1}n$$
$$2*z_p\sigma_T=z-w$$
$$\sigma_T=\frac{z-w}{2z_p}$$
We can use the quantile function of the normal distribution to find $\sigma_T$ as a function of $n$:
$$z_p=\sqrt2\text{erf}^{-1}(2p-1)$$
$$z_p=\sqrt2\text{erf}^{-1}\left(\frac{n-2}n\right)$$
$$\sigma_T=\frac{z-w}{2\sqrt2\text{erf}^{-1}\left(\frac{n-2}n\right)}$$
Or, we can just use quantile tables that give us $z_p$ for values of $p\approx p(n)$.
3. Bringing it all together
Our solution will be a particular Wiener Process distribution that behaves like the Brownian Bridge, but also adds up over time to the total distribution described in Part 2:
$$W(t)\sim B(t)$$
$$B_T=\int_{0}^{1}W(t)dt$$
A helpful property is that the product of two independent normal distributions is another normal distribution. So, we can decompose our desired model into the joint probability of $B(t)$ and an independent, time-invariant normal distribution $N_C$:
$$W(t)\ = N_CB(t)$$
$$B_T=N_C\int_{0}^{1}B(t)dt$$
The time-integral of a Brownian or Wiener process is another specific normal distribution:
$$B_T=N_CN(0,\tfrac13)$$
In general, the product of two normal distributions is the normal distribution with mean and variance as:
$$\mu=\frac{\mu_1\sigma^2_2+\mu_2\sigma^2_1}{\sigma^2_1+\sigma^2_2}$$
$$\sigma^2=\frac{\sigma^2_1\sigma^2_2}{\sigma^2_1+\sigma^2_2}$$
So then $N_C$ is fully specified by:
$$\sigma^2_T=\frac{\sigma^2_C\tfrac13}{\sigma^2_C+\tfrac13}$$
$$\sigma^2_C=\frac{\sigma^2_T}{1-3\sigma^2_T}$$
and
$$\mu_T=\frac{\tfrac13\mu_C+0}{\sigma^2_C+\tfrac13}$$
$$\mu_C=\mu_T(3\sigma^2_C+1)$$
$$\mu_C=\frac{\mu_T}{1-3\sigma^2_T}$$
Finally, we can use the same property to fully specify our model $W(t)=N_CB(t)=N(\mu_W,\sigma_W^2)$ in terms of the values solved for in Parts 1 and 2:
$$\mu_W=\frac{\mu_C\sigma^2_t+\mu_t\sigma^2_C}{\sigma^2_C+\sigma^2_t}$$
$$\sigma_W^2=\frac{\sigma^2_C\sigma^2_t}{\sigma^2_C+\sigma^2_t}$$
