R - Returning the full joint probability distribution of a dataset

Question

I have a data.frame with categorical variables as below:

bird.data <- data.frame( id =  rep(2,500),
                          colour = sample(c("Red", "Blue", "Yellow", "Green"), 500, replace=T, c(0.15,0.45, 0.20, 0.20)),
                          size = sample(c("Large", "Medium","Small"), 500, replace = T, c(0.33,0.33, 0.33)),
                          texture = sample(c("Hard", "Soft"), 500, replace = T, prob = c(0.55,0.45))
)

Is there an easy way to return the full joint distribution, P(colour,size,texture) of the dataset using R? For the dataset above, this would be a cube with dimensions: with(bird.data, levels(colour) * levels(size) * levels(texture)).

For instance, for the dataset given above, I want to be able to store every information as below within the cube:

# P(colour="Red", size="Small", texture= "Hard")
p_Red_Small_Hard <- nrow(bird.data[ bird.data$colour== "Red" & bird.data$size == "Small" & bird.data$texture =="Hard", ]) / nrow(bird.data)