Has anyone done the math on this and can share the result?
What is the average amount of rolls with different amount of dice and how does the distribution look?
It's been roughly 30 years since I last did this kind of calculations...
--------------------
Edit: See post #23 below for a much shorter way to do it.

Another thought, you could add a mechanic for characters to push their luck by willfully removing a die from the pool in return for some beneficial effect - a hint, advantage on a roll, rerolling damage, acquisition of a tangentially beneficial item, etc. Furthermore, you can also use gradiants...
www.enworld.org ---------------------I think the attached R code using a transition matrix is correct for giving you what turn you run out on.
#Set the number of starting dice
dat1<-6
#Set how many turns to check for
nturn<-50
#Set the number of sides on the dice
nside<-6
#Make the initial state vector
s1<-c(1,rep(0,dat1))
#Transition matrix
tmat<-matrix(0,ncol=dat1+1,nrow=dat1+1)
for (i in 1: (dat1+1)){
temp<-dbinom(((dat1-i+1):0),
dat1-i+1,(nside-1)/nside)
tmat[i,(i: (dat1+1))]<-temp
}
#Get probability it has run out by that turn
current<-s1%*%tmat
pendby<-current[dat1+1]
for (i in 1: (nturn-1)){
current<-current%*%tmat
pendby<-c(pendby,current[dat1+1])}
#Get probability it has run out on that turn
pendon<-pendby[1]
for (i in 2:nturn){
pendon<-c(pendon,pendby[ i] -pendby[i-1])}
pendon<-c(pendon,1-sum(pendon))
x<-cbind(1: (nturn+1),round(pendon,4),round(c(pendby,1),4))
colnames(x)=c("Turn","Prob On","Prob By")
x
The first value of pendon is the probability of running out of dice on turn 1, etc... The last value is the probability it takes more than whatever you set the max number of turns to check for. pendby will give you the probability of running out on that turn or before.
I checked it by writing one that did a simulation.
#Set the number of starting dice
dat1<-6
#Set how many turns to check for
nturn<-50
#Set the number of sides on the dice
nside<-6
#Set number of simulations
nsims<-100000
simcount<-rep(0,(nturn+1))
for (i in 1:nsims){
curcount<-dat1
turn<-0
while ((curcount>0)&(turn<(nturn+1))){
turn<-turn+1
curcount<-rbinom(1,curcount,(nside-1)/nside)
if (curcount==0){simcount[ i]<-turn}
if ((curcount>0)&(turn==nturn)){simcount[ i]<- (nturn+1)}
}}
table(simcount)/nsims
#Compare them to make sure it worked
cbind(1: (nturn+1),round(pendon,4),
round(table(simcount)/nsims,4))
Here are the probabilities for starting with three dice and going up to 20 turns. (Remember the last entry is "More than 20").
Turn Prob On Prob By
[1,] 1 0.0046 0.0046
[2,] 2 0.0239 0.0285
[3,] 3 0.0462 0.0748
[4,] 4 0.0640 0.1388
[5,] 5 0.0752 0.2140
[6,] 6 0.0802 0.2942
[7,] 7 0.0805 0.3747
[8,] 8 0.0773 0.4520
[9,] 9 0.0720 0.5240
[10,] 10 0.0655 0.5895
[11,] 11 0.0586 0.6481
[12,] 12 0.0517 0.6999
[13,] 13 0.0451 0.7450
[14,] 14 0.0391 0.7841
[15,] 15 0.0336 0.8176
[16,] 16 0.0287 0.8464
[17,] 17 0.0244 0.8708
[18,] 18 0.0207 0.8915
[19,] 19 0.0175 0.9090
[20,] 20 0.0148 0.9238
[21,] 21 0.0762 1.0000
Here it is for six dice
Turn Prob On Prob By
[1,] 1 0.0000 0.0000
[2,] 2 0.0008 0.0008
[3,] 3 0.0048 0.0056
[4,] 4 0.0137 0.0193
[5,] 5 0.0265 0.0458
[6,] 6 0.0408 0.0866
[7,] 7 0.0538 0.1404
[8,] 8 0.0639 0.2043
[9,] 9 0.0703 0.2746
[10,] 10 0.0730 0.3475
[11,] 11 0.0725 0.4201
[12,] 12 0.0697 0.4898
[13,] 13 0.0652 0.5550
[14,] 14 0.0597 0.6148
[15,] 15 0.0538 0.6685
[16,] 16 0.0478 0.7163
[17,] 17 0.0419 0.7583
[18,] 18 0.0365 0.7948
[19,] 19 0.0315 0.8263
[20,] 20 0.0271 0.8534
[21,] 21 0.1466 1.0000
Edit: I think I fixed all the auto-emoji and mark-up making in the code snippets.
ncG1vNJzZmivp6x7prrWqKmlnF6kv6h706GpnpmUqHyqutOrpp2tk567qHnToZxmm5%2Bqu7WwzrClZpyZmLJuucScn5qmmZh7d4WUamdrZw%3D%3D