Toss Up! is a very simple dice game that I’ve always wanted to simulate but never got around to doing so (until now!). This post outlines how to simulate a Toss Up! game in R, as well as how to evaluate the effectiveness of different game strategies. All the code for this blog post is available here.
Rules
The official rules for Toss Up! are available here. Here is an abbreviated version:
- As shown in the picture above, the game uses 10 six-sided dice. Each die has 3 green sides, 2 yellow sides and one red side.
- For each turn, a player rolls all 10 dice. If any greens are rolled, set them aside.
- At this point, the player can choose to end the turn and “bank” the points (one point for one green), or keep rolling the remaining dice.
- After every roll, greens are set aside. If you roll enough times to make all 10 dice green, you have 10 points and you can either bank or roll all 10 dice again.
- A player can keep going until he/she decides to stop rolling and bank all the points earned on this turn.
- Catch: If, on a roll, none of the dice are green and at least one is red, the turn ends and no points can be banked.
- First player to reach 100 points wins.
Simulating Toss Up!
There are several different ways to implement Toss Up!: I describe my (two-player) version below. To allow for greater flexibility for the resulting code, I implement the game with three global constants that the programmer can change:
NUM_DICE: The number of dice used in the game (default is 10).FACE_PROBS: A numeric vector of length 3 denoting the probability that a die comes up red, yellow or green respectively (default is c(1/6, 2/6, 3/6)).POINTS_TO_WIN: The number of points a player needs to obtain to win the game (default is 100).
Step 1: How can we describe the state of a game?
The state of the game can be encapsulated with 4 pieces of information:
current_player: who’s turn it is to make a decision.scores: a vector of length 2 containing the scores for players 1 and 2 respectively.turn_points: Number of points scored on the current turn so far (these points have not been banked yet).dice_rolls: A vector of variable length denoting the outcome of the dice rolls (0: red, 1: yellow, 2: green).
In my implementation, the state is stored as a list with the 4 elements above.
Step 2: Updating the state
From a given state, there are 3 possibilities for what comes next:
- The dice rolls don’t have any greens and at least one red. In this case, the current players turn is over. We need to change
current_player, reset turn_points to 0, and re-roll the dice (NUM_DICE of them).
- The dice rolls either (i) have at least one green, or (ii) have no reds. In this case, the current player has a choice of what to do.
- If the player chooses to bank, then we need to update
scores, reset turn_points to 0, change current_player and re-roll the dice (NUM_DICE of them).
- If the player chooses to roll, then we need to update
turn_points and re-roll just the dice that were not green.
The function updateState below does all of the above. I have also added a verbose option which, if set to TRUE, prints information on the game state to the console. (The function DiceToString is a small helper function for printing out the dice rolls.)
DiceToString <- function(dice_rolls) {
return(paste(sum(dice_rolls == 2), "Green,",
sum(dice_rolls == 1), "Yellow,",
sum(dice_rolls == 0), "Red"))
}
# Executes the current player's action and updates the state. Is also used if
# no greens and at least one red is rolled.
UpdateState <- function(state, action, verbose = FALSE) {
if (verbose) {
cat("Current roll:", DiceToString(state$dice_rolls))
if (action != "rolled at least 1 red and no green")
cat(" (bank", state$turn_points + sum(state$dice_rolls == 2), "pts?)",
fill = TRUE)
else
cat("", fill = TRUE)
cat(paste0("Player ", state$current_player, " ", action, "s"),
fill = TRUE)
}
if (action == "bank") {
# action = "bank": bank the points current player earned this turn, then
# re-roll the dice.
state$scores[state$current_player] <- state$scores[state$current_player] +
state$turn_points + sum(state$dice_rolls == 2)
state$turn_points <- 0
state$dice_rolls <- sample(0:2, size = NUM_DICE, replace = TRUE,
prob = FACE_PROBS)
state$current_player <- 3 - state$current_player
} else if (action == "roll") {
# action = "roll": add the green dice to points earned this turn, then
# re-roll the remaining dice.
state$turn_points <- state$turn_points + sum(state$dice_rolls == 2)
ndice <- sum(state$dice_rolls != 2)
if (ndice == 0) ndice <- NUM_DICE
state$dice_rolls <- sample(0:2, size = ndice, replace = TRUE,
prob = FACE_PROBS)
} else if (action == "rolled at least 1 red and no green") {
# action = "rolled at least 1 red and no green":
# no points banked, turn ends, re-roll dice.
state$turn_points <- 0
state$dice_rolls <- sample(0:2, size = NUM_DICE, replace = TRUE,
prob = FACE_PROBS)
state$current_player <- 3 - state$current_player
} else {
stop("action must be 'bank', 'roll', or 'rolled at least 1 red and no green'")
}
if (verbose) {
if (action != "roll") {
cat("Current scores:", state$scores, fill = TRUE)
cat("", fill = TRUE)
}
}
return(state)
}
Step 3: How to express player behavior and strategy?
We can think of a player as a black box which takes a game state as an input and outputs a decision, “bank” or “roll”. In other words, the player is a function!
Below are two extremely simple strategies expressed as functions. The first strategy simply banks after the first roll. The second strategy banks once more than target points can be earned from the turn.
# strategy that stops after one roll
OneRoll <- function(state) {
return("bank")
}
# strategy that stops only when the turn yields > `target` points
OverTargetPoints <- function(state, target = 10) {
if (state$turn_points + sum(state$dice_rolls == 2) > target) {
return("bank")
} else {
return("roll")
}
}
(Note: In my implementation, strategy functions should assume that they are current_player: that is how they can determine their current score and that of their opponent correctly.)
Step 4: Simulating a full game
We can now put the pieces from the previous steps together to simulate a full game of Toss Up!. The SimulateGame function takes in two strategy functions as input. It sets up the initial state, then updates the state repeatedly until the game ends.
SimulateGame <- function(Strategy1, Strategy2, verbose = FALSE) {
# set up initial state
state <- list(current_player = 1,
scores = c(0, 0),
turn_points = 0,
dice_rolls = sample(0:2, size = NUM_DICE, replace = TRUE,
prob = FACE_PROBS))
# check if no greens and at least one red, if so change player
while (sum(state$dice_rolls == 2) == 0 && sum(state$dice_rolls == 0) > 0) {
state <- UpdateState(state, "rolled at least 1 red and no green", verbose)
}
# while the game has not ended:
# - get the next action from the current player's strategy
# - update the state
while (max(state$scores) < POINTS_TO_WIN) {
if (state$current_player == 1) {
action <- Strategy1(state)
} else {
action <- Strategy2(state)
}
state <- UpdateState(state, action, verbose)
# check if no greens and at least one red
while (sum(state$dice_rolls == 2) == 0 && sum(state$dice_rolls == 0) > 0) {
state <- UpdateState(state, "rolled at least 1 red and no green", verbose)
}
}
# game has ended: return winner
if (verbose) {
cat(paste("Game ends: Player", which.max(state$scores), "wins!"),
fill = TRUE)
}
return(which.max(state$scores))
}
Two examples of simulated games
The code below shows what a simulated game of Toss Up! might look like. In the code snippet below, players need 20 points to win. The first player uses the super conservative strategy of banking immediately, while the second player tries to win the entire game in one turn.
NUM_DICE <- 10
FACE_PROBS <- c(1/6, 2/6, 3/6)
POINTS_TO_WIN <- 20
set.seed(1)
winner <- SimulateGame(OneRoll,
function(state) { OverTargetPoints(state, 19) },
verbose = TRUE)
# Current roll: 4 Green, 3 Yellow, 3 Red (bank 4 pts?)
# Player 1 banks
# Current scores: 4 0
#
# Current roll: 5 Green, 4 Yellow, 1 Red (bank 5 pts?)
# Player 2 rolls
# Current roll: 3 Green, 1 Yellow, 1 Red (bank 8 pts?)
# Player 2 rolls
# Current roll: 2 Green, 0 Yellow, 0 Red (bank 10 pts?)
# Player 2 rolls
# Current roll: 5 Green, 4 Yellow, 1 Red (bank 15 pts?)
# Player 2 rolls
# Current roll: 2 Green, 3 Yellow, 0 Red (bank 17 pts?)
# Player 2 rolls
# Current roll: 0 Green, 3 Yellow, 0 Red (bank 17 pts?)
# Player 2 rolls
# Current roll: 2 Green, 1 Yellow, 0 Red (bank 19 pts?)
# Player 2 rolls
# Current roll: 0 Green, 1 Yellow, 0 Red (bank 19 pts?)
# Player 2 rolls
# Current roll: 0 Green, 1 Yellow, 0 Red (bank 19 pts?)
# Player 2 rolls
# Current roll: 1 Green, 0 Yellow, 0 Red (bank 20 pts?)
# Player 2 banks
# Current scores: 4 20
#
# Game ends: Player 2 wins!
In this particular simulation, it paid off to be a bit more risk taking. As the simulation below shows, that’s not always the case.
NUM_DICE <- 10
FACE_PROBS <- c(0.5, 0.0, 0.5)
POINTS_TO_WIN <- 20
set.seed(1)
winner <- SimulateGame(OneRoll,
function(state) { OverTargetPoints(state, 19) },
verbose = TRUE)
# Current roll: 6 Green, 0 Yellow, 4 Red (bank 6 pts?)
# Player 1 banks
# Current scores: 6 0
#
# Current roll: 5 Green, 0 Yellow, 5 Red (bank 5 pts?)
# Player 2 rolls
# Current roll: 2 Green, 0 Yellow, 3 Red (bank 7 pts?)
# Player 2 rolls
# Current roll: 0 Green, 0 Yellow, 3 Red
# Player 2 rolled at least 1 red and no greens
# Current scores: 6 0
#
# Current roll: 5 Green, 0 Yellow, 5 Red (bank 5 pts?)
# Player 1 banks
# Current scores: 11 0
#
# Current roll: 7 Green, 0 Yellow, 3 Red (bank 7 pts?)
# Player 2 rolls
# Current roll: 2 Green, 0 Yellow, 1 Red (bank 9 pts?)
# Player 2 rolls
# Current roll: 1 Green, 0 Yellow, 0 Red (bank 10 pts?)
# Player 2 rolls
# Current roll: 3 Green, 0 Yellow, 7 Red (bank 13 pts?)
# Player 2 rolls
# Current roll: 2 Green, 0 Yellow, 5 Red (bank 15 pts?)
# Player 2 rolls
# Current roll: 2 Green, 0 Yellow, 3 Red (bank 17 pts?)
# Player 2 rolls
# Current roll: 2 Green, 0 Yellow, 1 Red (bank 19 pts?)
# Player 2 rolls
# Current roll: 0 Green, 0 Yellow, 1 Red
# Player 2 rolled at least 1 red and no greens
# Current scores: 11 0
#
# Current roll: 5 Green, 0 Yellow, 5 Red (bank 5 pts?)
# Player 1 banks
# Current scores: 16 0
#
# Current roll: 4 Green, 0 Yellow, 6 Red (bank 4 pts?)
# Player 2 rolls
# Current roll: 4 Green, 0 Yellow, 2 Red (bank 8 pts?)
# Player 2 rolls
# Current roll: 1 Green, 0 Yellow, 1 Red (bank 9 pts?)
# Player 2 rolls
# Current roll: 0 Green, 0 Yellow, 1 Red
# Player 2 rolled at least 1 red and no greens
# Current scores: 16 0
#
# Current roll: 5 Green, 0 Yellow, 5 Red (bank 5 pts?)
# Player 1 banks
# Current scores: 21 0
#
# Game ends: Player 1 wins!
Comparing some simple strategies
We can’t tell whether one strategy is better than another by looking at a single game since there is so much randomness involved. What we should do is simulate many games and see which strategy wins out over the long run.
OneRoll vs. OneRoll
Let’s do a sanity check. Here, we run 10,000 games of the strategy OneRoll vs. itself (this takes around 4 seconds on my machine). Since the strategy for both players is the same, we might expect player 1 to win around 50% of the time right?
NUM_DICE <- 10
FACE_PROBS <- c(1/6, 2/6, 3/6)
POINTS_TO_WIN <- 100
nsim <- 10000
set.seed(1)
winners <- replicate(nsim, SimulateGame(OneRoll, OneRoll))
table(winners) # not 50-50: player who starts first has an advantage
# winners
# 1 2
# 5910 4090
It looks like player 1 wins around 59% of the time! On second thought, we would expect player 1 to win more because he/she has a first mover advantage: if player 1 reaches 100 points, he/she wins even if player 2 might reach 100 points on the very next turn.
To make this sanity check work, we should have the players each go first half the time. The code below achieves that:
# for even numbered simulations, let player 2 start first
winners2 <- winners
winners2[2 * 1:(nsim/2)] <- 3 - winners2[2 * 1:(nsim/2)]
table(winners2)
# winners2
# 1 2
# 5030 4970
Now player 1 wins just 50.3% of the time, which is much closer to the expected 50%.
OneRoll vs. >20 points
Next, let’s compare the “one roll” strategy with the strategy which stops once the player can bank more than 20 points for the turn:
set.seed(1)
winners <- replicate(nsim, SimulateGame(
OneRoll,
function(state) { OverTargetPoints(state, 20) }))
table(winners)
# winners
# 1 2
# 75 9925
Wow! The strategy of banking only when >20 points have been earned on a turn wins almost all the time, even though it doesn’t have the first mover advantage!
>20 points vs. >50 points
Taking a bit more risk than always banking seems to pay off. What about even more risk? How does banking after >20 points compare with banking after >50 points?
set.seed(1)
winners <- replicate(nsim, SimulateGame(
function(state) { OverTargetPoints(state, 20) },
function(state) { OverTargetPoints(state, 50) }))
table(winners)
# winners
# 1 2
# 6167 3833
# switch order of players
set.seed(1)
winners <- replicate(nsim, SimulateGame(
function(state) { OverTargetPoints(state, 50) },
function(state) { OverTargetPoints(state, 20) }))
table(winners)
# winners
# 1 2
# 4414 5586
The >20 points strategy won 61.7% of the games when it went first, and won 55.9% of the games when it went second, so it’s clearly a superior strategy.
Where can we go from here?
There are several future directions we could take with this code base.
- The strategies presented here are very simple. Can we program more complex strategies (e.g. some that depend on the score of the other player, or on how close the player is to winning)?
- Among strategies of the form “bank once you have more than xx points”, is there an optimal value for xx?
- What is the optimal strategy for Toss Up!? Can we learn it through methods such as Q-learning?
- Is there even such a thing as an optimal strategy? Do there exist 3 strategies such that A beats B, B beats C, but C beats A?
- Can we quantify the first mover advantage?
- Can we modify the code so that we can visualize the players’ scores as the game progresses?
- Can we modify the code so that more than 2 players can play?
, there are constants
and
such that
are independent copies of
, a sequence of positive numbers
and a sequence of numbers
such that
denotes convergence in distribution.)
,
,
,
such that its characteristic function has the following form:
where
, we get the normal distribution with mean
and variance
.
and
, we get the Cauchy distribution with scale parameter
and shift parameter
and
, we get the Lévy distribution with scale parameter
with non-negative entries is said to be doubly stochastic if each of its rows and columns sum to 1, i.e.
are both doubly stochastic, then their product 
is also doubly stochastic. We can prove this directly by some algebra. For 
‘s rows: for each 
,
,
and 
are both doubly stochastic, we can view them as transition probability matrices for Markov chains denoted 
and 
. If we think of a new process where each step consists of taking one move according to 
is right stochastic, implying that 
is said to mean independent of ![\mathbb{E}[Y \mid X = x] = \mathbb{E}[Y]](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+x%5D+%3D+%5Cmathbb%7BE%7D%5BY%5D&bg=ffffff&fg=333333&s=0&c=20201002)
such that the probability mass/density of ![\mathbb{E}[Y \mid X]](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+%5Cmid+X%5D&bg=ffffff&fg=333333&s=0&c=20201002)
is the unique random variable such that for any 
(the 
-algebra generated by ![\mathbb{E}[\mathbb{E}[Y \mid X] 1_A] = \mathbb{E}[Y 1_A].](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5B%5Cmathbb%7BE%7D%5BY+%5Cmid+X%5D+1_A%5D+%3D+%5Cmathbb%7BE%7D%5BY+1_A%5D.&bg=ffffff&fg=333333&s=0&c=20201002)
![\mathbb{E}[Y 1_A] = \mathbb{E}[Y] \mathbb{E}[1_A] = \mathbb{E}[\mathbb{E}[Y]1_A].](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+1_A%5D+%3D+%5Cmathbb%7BE%7D%5BY%5D+%5Cmathbb%7BE%7D%5B1_A%5D+%3D+%5Cmathbb%7BE%7D%5B%5Cmathbb%7BE%7D%5BY%5D1_A%5D.&bg=ffffff&fg=333333&s=0&c=20201002)
![E[Y]](https://s0.wp.com/latex.php?latex=E%5BY%5D&bg=ffffff&fg=333333&s=0&c=20201002)
, i.e. 
is the same for all ![\mathbb{E}[Y \mid X = x]](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+x%5D&bg=ffffff&fg=333333&s=0&c=20201002)
is a function of the conditional distribution 
, it follows that ![\mathbb{E}[Y]](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY%5D&bg=ffffff&fg=333333&s=0&c=20201002)
.)![\begin{aligned} \mathbb{E}[XY] &= \mathbb{E}[\mathbb{E}[XY \mid X]] \\ &= \mathbb{E}[X \mathbb{E}[Y \mid X]] \\ &= \mathbb{E}[X \mathbb{E}[Y]] \\ &= \mathbb{E}[Y] \mathbb{E}[X], \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+%5Cmathbb%7BE%7D%5BXY%5D+%26%3D+%5Cmathbb%7BE%7D%5B%5Cmathbb%7BE%7D%5BXY+%5Cmid+X%5D%5D+%5C%5C++%26%3D+%5Cmathbb%7BE%7D%5BX+%5Cmathbb%7BE%7D%5BY+%5Cmid+X%5D%5D+%5C%5C++%26%3D+%5Cmathbb%7BE%7D%5BX+%5Cmathbb%7BE%7D%5BY%5D%5D+%5C%5C++%26%3D+%5Cmathbb%7BE%7D%5BY%5D+%5Cmathbb%7BE%7D%5BX%5D%2C++%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
.
and 
, with ![\theta \sim \text{Unif}[0, 2\pi]](https://s0.wp.com/latex.php?latex=%5Ctheta+%5Csim+%5Ctext%7BUnif%7D%5B0%2C+2%5Cpi%5D&bg=ffffff&fg=333333&s=0&c=20201002)
. For any ![x \in [-1, 1]](https://s0.wp.com/latex.php?latex=x+%5Cin+%5B-1%2C+1%5D&bg=ffffff&fg=333333&s=0&c=20201002)
, 
, so ![\mathbb{E}[Y \mid X = x] = 0](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+x%5D+%3D+0&bg=ffffff&fg=333333&s=0&c=20201002)
, i.e. 
and 
, but since 
,
with equal probability. With this setup,![\mathbb{E}[XY] = \mathbb{E}[X] = \mathbb{E}[Y] = 0](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BXY%5D+%3D+%5Cmathbb%7BE%7D%5BX%5D+%3D+%5Cmathbb%7BE%7D%5BY%5D+%3D+0&bg=ffffff&fg=333333&s=0&c=20201002)
,![\mathbb{E}[Y \mid X = 1] = 3, \quad \mathbb{E}[Y \mid X = -3] = 1,](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+1%5D+%3D+3%2C+%5Cquad+%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+-3%5D+%3D+1%2C&bg=ffffff&fg=333333&s=0&c=20201002)
with equal probability. With this setup,![\mathbb{E}[X \mid Y = 1] = \mathbb{E}[X \mid Y = -1] = 0,](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BX+%5Cmid+Y+%3D+1%5D+%3D+%5Cmathbb%7BE%7D%5BX+%5Cmid+Y+%3D+-1%5D+%3D+0%2C&bg=ffffff&fg=333333&s=0&c=20201002)
![\mathbb{E}[Y \mid X = 2] = 1, \quad \mathbb{E}[Y \mid X = 5] = -1,](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+2%5D+%3D+1%2C+%5Cquad+%5Cmathbb%7BE%7D%5BY+%5Cmid+X+%3D+5%5D+%3D+-1%2C&bg=ffffff&fg=333333&s=0&c=20201002)
, let ![U_{n,a} = \mathbb{E}\left[\dfrac{1}{a+X} \right]](https://s0.wp.com/latex.php?latex=U_%7Bn%2Ca%7D+%3D+%5Cmathbb%7BE%7D%5Cleft%5B%5Cdfrac%7B1%7D%7Ba%2BX%7D+%5Cright%5D&bg=ffffff&fg=333333&s=0&c=20201002)
with 
. Then
![\mathbb{E}\left[\dfrac{1}{a+X} \right]](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5Cleft%5B%5Cdfrac%7B1%7D%7Ba%2BX%7D+%5Cright%5D&bg=ffffff&fg=333333&s=0&c=20201002)
, we can write out an expression for 
by iteratively substituting the recurrence relation into itself. (Maybe we can use that to get a single explicit formula; I didn’t try too hard.)
,![\begin{aligned} \mathbb{E}\left[ \dfrac{1}{a + X} \right] = \int_0^1 t^{a-1} \cdot P_X(t) dt, \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+%5Cmathbb%7BE%7D%5Cleft%5B+%5Cdfrac%7B1%7D%7Ba+%2B+X%7D+%5Cright%5D+%3D+%5Cint_0%5E1+t%5E%7Ba-1%7D+%5Ccdot+P_X%28t%29+dt%2C+%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
is the 
, we have 
for 
,![\begin{aligned} U_{n, 1} = \int_0^1 t^0 (q + pt)^n dt = \left[ \dfrac{(q+pt)^{n+1}}{(n+1)p} \right]_0^1 = \dfrac{1 - q^{n+1}}{(n+1)p}. \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+U_%7Bn%2C+1%7D+%3D+%5Cint_0%5E1+t%5E0+%28q+%2B+pt%29%5En+dt++%3D+%5Cleft%5B+%5Cdfrac%7B%28q%2Bpt%29%5E%7Bn%2B1%7D%7D%7B%28n%2B1%29p%7D+%5Cright%5D_0%5E1++%3D+%5Cdfrac%7B1+-+q%5E%7Bn%2B1%7D%7D%7B%28n%2B1%29p%7D.+%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
, using integration by parts,![\begin{aligned} U_{n, a} &= \int_0^1 t^{a-1} (q + pt)^n dt \\ &= \left[ t^{a-1} \cdot \dfrac{(q+pt)^{n+1}}{(n+1)p} \right]_0^1 - \int_0^1 \dfrac{(q+pt)^{n+1}}{(n+1)p} \cdot (a-1)t^{a-2} dt \\ &= \left[ \dfrac{1}{(n+1)p} - 0 \right] - \dfrac{a-1}{(n+1)p} \int_0^1 t^{(a-1)-1} (q+pt)^{n+1} dt \\ &= \dfrac{1 - (a-1) U_{n+1, a-1}}{(n+1)p}, \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+U_%7Bn%2C+a%7D+%26%3D+%5Cint_0%5E1+t%5E%7Ba-1%7D+%28q+%2B+pt%29%5En+dt+%5C%5C++%26%3D+%5Cleft%5B+t%5E%7Ba-1%7D+%5Ccdot+%5Cdfrac%7B%28q%2Bpt%29%5E%7Bn%2B1%7D%7D%7B%28n%2B1%29p%7D+%5Cright%5D_0%5E1+-+%5Cint_0%5E1+%5Cdfrac%7B%28q%2Bpt%29%5E%7Bn%2B1%7D%7D%7B%28n%2B1%29p%7D+%5Ccdot+%28a-1%29t%5E%7Ba-2%7D+dt+%5C%5C++%26%3D+%5Cleft%5B+%5Cdfrac%7B1%7D%7B%28n%2B1%29p%7D+-+0+%5Cright%5D+-+%5Cdfrac%7Ba-1%7D%7B%28n%2B1%29p%7D+%5Cint_0%5E1+t%5E%7B%28a-1%29-1%7D+%28q%2Bpt%29%5E%7Bn%2B1%7D+dt+%5C%5C++%26%3D+%5Cdfrac%7B1+-+%28a-1%29+U_%7Bn%2B1%2C+a-1%7D%7D%7B%28n%2B1%29p%7D%2C++%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
exactly with the recurrence relation, while the second does a Monte Carlo simulation. Here is the code that verifies the identity for 
, 
and 
:
, and let ![\begin{aligned} P_X(t) = \mathbb{E}[t^X] = \sum_{x=0}^\infty t^x \mathbb{P} \{ X = x\}. \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+P_X%28t%29+%3D+%5Cmathbb%7BE%7D%5Bt%5EX%5D+%3D+%5Csum_%7Bx%3D0%7D%5E%5Cinfty+t%5Ex+%5Cmathbb%7BP%7D+%5C%7B+X+%3D+x%5C%7D.+%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
![\begin{aligned} \mathbb{E}\left[ \dfrac{1}{a + X} \right] = \int_0^1 t^{a-1} \cdot P_X(t) dt. \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+%5Cmathbb%7BE%7D%5Cleft%5B+%5Cdfrac%7B1%7D%7Ba+%2B+X%7D+%5Cright%5D+%3D+%5Cint_0%5E1+t%5E%7Ba-1%7D+%5Ccdot+P_X%28t%29+dt.+%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
![\mathbb{E}\left[ \frac{1}{1+X} \right]](https://s0.wp.com/latex.php?latex=%5Cmathbb%7BE%7D%5Cleft%5B+%5Cfrac%7B1%7D%7B1%2BX%7D+%5Cright%5D&bg=ffffff&fg=333333&s=0&c=20201002)
for ![\begin{aligned} \int_0^1 t^{a-1} \cdot P_X(t) dt &= \int_0^1 t^{a-1} \cdot \left( \sum_{x=0}^\infty t^x \mathbb{P} \{ X = x\} \right) dt \\ &= \sum_{x=0}^\infty \mathbb{P} \{ X = x\} \left( \int_0^1 t^{x + a - 1} dt \right) \\ &= \sum_{x=0}^\infty \mathbb{P} \{ X = x\} \left[ \dfrac{t^{x+a}}{x+a} \right]_0^1 \\ &= \sum_{x=0}^\infty \mathbb{P} \{ X = x\} \dfrac{1}{x+a} \\ &= \mathbb{E}\left[\dfrac{1}{X+a} \right]. \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+%5Cint_0%5E1+t%5E%7Ba-1%7D+%5Ccdot+P_X%28t%29+dt+%26%3D+%5Cint_0%5E1+t%5E%7Ba-1%7D+%5Ccdot+%5Cleft%28+%5Csum_%7Bx%3D0%7D%5E%5Cinfty+t%5Ex+%5Cmathbb%7BP%7D+%5C%7B+X+%3D+x%5C%7D+%5Cright%29+dt+%5C%5C++%26%3D+%5Csum_%7Bx%3D0%7D%5E%5Cinfty+%5Cmathbb%7BP%7D+%5C%7B+X+%3D+x%5C%7D+%5Cleft%28+%5Cint_0%5E1+t%5E%7Bx+%2B+a+-+1%7D+dt+%5Cright%29+%5C%5C++%26%3D+%5Csum_%7Bx%3D0%7D%5E%5Cinfty+%5Cmathbb%7BP%7D+%5C%7B+X+%3D+x%5C%7D+%5Cleft%5B+%5Cdfrac%7Bt%5E%7Bx%2Ba%7D%7D%7Bx%2Ba%7D+%5Cright%5D_0%5E1+%5C%5C++%26%3D+%5Csum_%7Bx%3D0%7D%5E%5Cinfty+%5Cmathbb%7BP%7D+%5C%7B+X+%3D+x%5C%7D+%5Cdfrac%7B1%7D%7Bx%2Ba%7D+%5C%5C++%26%3D+%5Cmathbb%7BE%7D%5Cleft%5B%5Cdfrac%7B1%7D%7BX%2Ba%7D+%5Cright%5D.+%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
distribution is through the question: “What is the probability that the random variable is at least x standard deviations (SDs) away from the mean?” For the most part, the smaller the degrees of freedom, the larger this probability was (more “heavy-tailed”), until we realized that the trend reversed for really small degrees of freedom (2.1 in the post). In fact, for 
, the variance of the 
with 
against 
? If I did that, I would have ended up with the plot below (for the log of the probabilities):
distribution is more heavy-tailed than the 
distribution??
. For each random variable 
, find the scale factor 
. For this value of 
and 
with 
, or a 1-in-1,000 chance.
: it seems that for any number of SDs, the probability of being that number of SD out for 
. Does that mean that 
degrees of freedom has SD 
. For 
, the SD is about 1.73 while for 
the SD is about 4.58, much larger! Taking this to the extreme, consider a 
. The variance is infinite in this case, so the random variable always takes values within 1 SD of the mean! Does it mean that this distribution is less heavy tailed than the normal distribution?
and parameters 
, commonly written as a vector 
. 
has the 
) if it has support on the probability simplex 
and its probability density function (PDF) satisfies
distribution is special because its PDF is constant over its support:
![U_1, \dots, U_{K-1} \sim \text{Unif}[0,1]](https://s0.wp.com/latex.php?latex=U_1%2C+%5Cdots%2C+U_%7BK-1%7D+%5Csim+%5Ctext%7BUnif%7D%5B0%2C1%5D&bg=ffffff&fg=333333&s=0&c=20201002)
, and let 
denote the order statistics for 
. Then
.)
. Then the joint density of the order statistics 
is
, we have
is the 
when 
.
, we have
until just two elements are left.
. If the random variables with subscripts 
and 
. Compute the PDF of 
directly:
,![\begin{aligned} f(x_1', x_3, \dots, x_K) &= \prod_{i=3}^K x_i^{\alpha_i - 1} \int_0^1 (x_1' v)^{\alpha_1 - 1} [x_1'(1 - v)]^{\alpha_2 - 1} (x_1')dv \\ &= (x_1')^{\alpha_1 + \alpha_2 - 1} \prod_{i=3}^K x_i^{\alpha_i - 1} \int_0^1 v^{\alpha_1 - 1} (1 - v)^{\alpha_2 - 1} dv \\ &\propto (x_1')^{\alpha_1 + \alpha_2 - 1} \prod_{i=3}^K x_i^{\alpha_i - 1}, \end{aligned}](https://s0.wp.com/latex.php?latex=%5Cbegin%7Baligned%7D+f%28x_1%27%2C+x_3%2C+%5Cdots%2C+x_K%29+%26%3D+%5Cprod_%7Bi%3D3%7D%5EK+x_i%5E%7B%5Calpha_i+-+1%7D+%5Cint_0%5E1+%28x_1%27+v%29%5E%7B%5Calpha_1+-+1%7D+%5Bx_1%27%281+-+v%29%5D%5E%7B%5Calpha_2+-+1%7D+%28x_1%27%29dv+%5C%5C++%26%3D+%28x_1%27%29%5E%7B%5Calpha_1+%2B+%5Calpha_2+-+1%7D+%5Cprod_%7Bi%3D3%7D%5EK+x_i%5E%7B%5Calpha_i+-+1%7D+%5Cint_0%5E1+v%5E%7B%5Calpha_1+-+1%7D+%281+-+v%29%5E%7B%5Calpha_2+-+1%7D+dv+%5C%5C++%26%5Cpropto+%28x_1%27%29%5E%7B%5Calpha_1+%2B+%5Calpha_2+-+1%7D+%5Cprod_%7Bi%3D3%7D%5EK+x_i%5E%7B%5Calpha_i+-+1%7D%2C+%5Cend%7Baligned%7D&bg=ffffff&fg=333333&s=0&c=20201002)
distribution.
Statistical Odds & Ends 