R, 198 bytes
function(m){r=c();a=diag(m);a[,1]=1;for(i in 3:m)for(j in 2:(i-1))a[i,j]=a[i-1,j-1]+a[i-j,j];p=sample(sum(a[m,]),1);while(m>0){b=a[m,];c=cumsum(b);x=min(which(c>=p));r=c(r,x);p=p-c[x]+b[x];m=m-x};r}
Ungolfed:
f <- function(m) {
r <- c()
a <- diag(m)
a[, 1] <- 1
for (i in 3:m)
for (j in 2:(i-1))
a[i, j] <- a[i-1, j-1] + a[i-j, j]
p <- sample(sum(a[m, ]), 1)
while (m > 0) {
b <- a[m, ]
c <- cumsum(b)
x <- min(which(c >= p))
r <- c(r, x)
p <- p - c[x] + b[x]
m <- m - x
}
return(r)
}
It follows the same structure as @dcsohl's great solution in Octavegreat solution in Octave, and is thus also based on the Knuth extract posted by @feersum.
I'll edit this later if I can come up with a more creative solution in R. In the meantime, any input is of course welcome.