R, 91 85 84 bytes

Edit: -1 byte thanks to pajonk

function(v,l=sum(v|T),`?`=sample,m=(l:0)[!(0:(2^l-2)?1)%/%2^(0:l)%%2]){v[m]=v[?m];v}

Try it online!

TIO link includes (reciprocal) distribution of shuffle_subsequence(1:3) to check conformity to tsh's calculations.

R, 91 85 84 82 bytes

Edit: -1 byte thanks to pajonk, then -2 bytes thanks to Giuseppe

function(v,l=sum(v|T),`?`=sample,m=(l:0)[!(2:2^l-2?1)%/%2^(0:l)%%2]){v[m]=v[?m];v}

Try it online!

TIO link includes (reciprocal) distribution of shuffle_subsequence(1:3) to check conformity to tsh's calculations.

R, 91 85 bytes

function(v,l=sum(v|T),s=sample,m=(l:0)[!s(0:(2^l-2),1)%/%2^(0:l)%%2]){v[m]=v[s(m)];v}

Try it online!

TIO link includes (reciprocal) distribution of shuffle_subsequence(1:3) to check conformity to tsh's calculations.

R, 91 85 84 bytes

Edit: -1 byte thanks to pajonk

function(v,l=sum(v|T),`?`=sample,m=(l:0)[!(0:(2^l-2)?1)%/%2^(0:l)%%2]){v[m]=v[?m];v}

Try it online!

TIO link includes (reciprocal) distribution of shuffle_subsequence(1:3) to check conformity to tsh's calculations.

R, 91 bytes

function(v,l=sum(v|T),s=sample){if(sum(m<-!s(0:(2^l-2),1)%/%2^(1:l-1)%%2)>1)v[m]=s(v[m]);v}

Try it online!

TIO link includes (reciprocal) distribution of shuffle_subsequence(1:3) to check conformity to tsh's calculations.

R, 91 85 bytes

function(v,l=sum(v|T),s=sample,m=(l:0)[!s(0:(2^l-2),1)%/%2^(0:l)%%2]){v[m]=v[s(m)];v}

Try it online!

TIO link includes (reciprocal) distribution of shuffle_subsequence(1:3) to check conformity to tsh's calculations.