A humble aspiring junior dev's guide to probabilistic thought. A series of rants and odd thoughts 2/4

What are the odds? I've spent quite a few hours invested in implementing a seemingly innocent feature of a game I'm currently working on as part of my learning in the world of functional programming with Elixir. I won't bore you with the details, I haven't even finished making the game, the prelude to discovering this problem is an innocent one. It's a terminal dungeon crawling game, where you venture into dark rooms in search of treasure and an eventual (hopeful) exit.

If you know anything about me? I loooovvveee games and RPGs especially JRPGs, I love-hate dungeon crawling and puzzles so much so that I write about it [NSFW!] so I've been having fun making this and curious about how all the gears fit, it is especially interesting to note how this relates to the modelling and simulation of objective functions something I'm quite interested in!!! Let's dive in and explore!

The probability of choosing a room is controlled by a call to Enum.random(iterable) which selects a random element from the iterable, using Erlang's internal pseudo-random generation algorithm. At this point in time, only three rooms existed (a room leading to the exit, A room full of monsters and another room full of monsters.)
Therefore it's safe to assume the mathematical probability of finding any room is P(1/3). This isn't so great is it? You could just start the game and one-third of the time, without experiencing even a single battle :( Oh I know the author suggests... let's add different probabilities to the rooms, such that rooms with battles appear more likely and rooms with an exit less likely. Okay... so I tried it, after the first hour or so I came up with a satisfactory(to me I thought) solution. This is it. I first added a probability key, with a value of an array of atoms to the struct.

defmodule DungeonCrawl.Room do  @moduledoc """  data structure representing a "room"  a room has "actions" borrowed  """  alias DungeonCrawl.Room, as: Room  alias DungeonCrawl.Room.Triggers, as: Triggers  import DungeonCrawl.Room.Action  defstruct name: nil, description: nil, actions: [], trigger: nil, probability: nil  def all,    do: [      %Room{        name: :exit,        description: "You can see a small light from a crack in the walls",        actions: [forward()],        trigger: Triggers.Exit,        probability: [:exit]      },    %Room{        name: :goblin,        description: "You can see an enemy blocking your path",        actions: [forward()],        trigger: Triggers.Enemy,        probability: [:goblin, :goblin, :goblin, :goblin, :goblin]      },    %Room{        name: :ogre,        description: "Something moves around in the dark, what do you do?",        actions: [forward()],        trigger: Triggers.Enemy,        probability: [:ogre, :ogre, :ogre]      },    ]end

At once you can probably see the naiveness of this solution. What computed the algorithm is this

 defp bias_probability(rooms) do    rooms    |> Enum.map(fn room -> room.probability end)    |> Enum.reduce(fn room, acc -> room ++ acc end)    |> Enum.random()  end

and finally now that I know the biased random outcome-

    rooms    |> Enum.find(fn room -> room.name == bias_probability(rooms) end)

Technically? This is a mathematically correct solution. It will compute a random value with a biased probability of P(0.1, 0.5, 0.3). I felt happy with myself, turned off my laptop and went to sleep! That was that! Or so I thought.

See there are two fundamental flaws of this approach:

1. I'm hard coding the probabilities, although modifying the probability of a room is possible, it will result in boiler plate code and unnecessary list traversal

2. The algorithm is not performant(in space complexity) atoms are not garbage collected in Elixir, and it would be unwise to arbitrarily create them on a whim! suppose I had a thousand rooms in the future for example

Sure enough, extra feature requirements came, the addition of a "difficulty level" that will bias this probability dynamically. So.. began my search, with a renewed courage to try mathematical problems I tweaked and tweaked, no luck. I asked my friends in engineering from school no luck! and then sure enough, a friend of mine with superior googling skill found a long lost blog post from 2017. This is what really helped me, without it I was lost in meaningless scribbles of probability theory and conditional problems. If you're interested in learning more about implementing a weighted probability algorithm please check out David Hills's blog it's really detailed and well done! I'll focus on implementing two algorithmic approaches I learned and adapting it to the context of this game. The linear search and one of the alias algorithms.

Probabilities can be unfair!

Linear scan the O(kay) method

O(n) runtime

A lovely idea is to search the linked list for an index value based on probable outcomes, wish I thought of this!! This will work by generating a random number between 0 - 1 and multiply that by the total probability distribution(For my use case, I'm assuming the list is sorted), then traverse the probability distribution subtracting each probability from the random probability, if it goes below zero? We've exhausted the probability distribution and found the index we need. Lovely solution. This works because we mathematically assume random.random() will generate a perfectly random number(spoiler it doesn't, but we don't care). If the probability is really low? It will hit the lowest index less often, and if the probability is high? It will hit that index more often. Run the code below to get a sense of it if you like.

import randomdef weighted_random(weights):    """    INPUT probabilities - [0.2, 0.3, 0.5]    OUTPUT index of random_weight - [0, 1, 2]    """    remaining_distance = random.random() * sum(weights)    # probability distribution sample size * random integer    for (i, weight) in enumerate(weights):        # [{0, 0.2}, {1, 0.3}, {2, 0.5}]        remaining_distance -= weights[i]        #print("debug", remaining_distance)        if remaining_distance < 0:            return i# Repeat the experiment to observe bias factorfor e in range(10):  print(f"exp trial{e}",weighted_random([0.2, 0.3, 0.5]))

This seems really cool doesn't it? I can represent the P(n) of my game as floats and dynamically update them, keeping the purity of the function. A quick and dirty Elixir version of this looks like

defmodule Prob do  def bias_probability(weights) do    # [i[0] P=0.2, i[1] P=0.5 , i[2] P=0.3] s.s = 3    # P(n) = distribution sample len bias * random() - P[i] < 0    len = Enum.count(weights)    # could increase range for greater float precision as you like    distance = Enum.random(1..99)/100 * Enum.sum(weights)    # bias factor    enumerated_list = Enum.zip(0..len, weights)    Enum.reduce_while(enumerated_list, distance, fn (_weight = {i, w}, acc) ->      if (acc - w) > 0, do: {:cont, acc - w}, else: {:halt, i}    end)  endendIO.inspect Prob.bias_probability([0.2, 0.3, 0.5])

and then to find the biased index.

    rooms    |> Enum.map(fn room -> room.probability end)    |> Enum.fetch(bias_probability(rooms))

but now that I've gone through the trouble of investigating weighted probability, I don't wanna settle for a linear search. Might as well make the game with the fastest possible approach. Let's go back in history, and take a look at one of the alias algorithms.

Aliasing(Vose) the O yes! method

O(n) alias table prep + O(1) lookup

This one is a bit tricky to understand at first. If you're interested in a much more technical deep dive into the various approaches for solving this problem, please check out this article by Keith Schwarz it's dense, comprehensive and technical. Onto aliasing!

The idea is honestly very clever! and in my opinion counter-intuitive, which is what makes it's constant-time runtime lookup so impressive. Let's take a step back, before diving in. An alias is like a nickname we give to stuff right? You have a friend called Samuel, and you call him Sam. His name isn't Sam and he hates that nickname, but alas, it stuck... poor Samuel, but if you see his face anywhere you know Sam is Samuel because you're such good friends. The alias method is a lot like this in principle. You have a biased probability weight distribution and it approaches some limit and that limit can be used to generate a custom distribution. Let's talk about Sam for a bit.

See Sam's nickname is called an alias table and your friendship(how you know Sam is Samuel) is the algorithm. Now that we have some hint as to what we're doing, here's the python implementation with my comments from Bruce Hill's blog. Scan through it, maybe run it if you're adventurous. We'll break it down step by step next. (If you're familiar with Java check out Keith Schwarz's implementation instead

def prepare_aliased_randomizer(weights):    N = len(weights) # limit    avg = sum(weights)/N # alias partition    aliases = [(1, None)] * N    # [(1, None), (1, None), (1, None)] --> alias table    # Bucketting (pseudo quick sort like behaviour*)    # weight/avg < | > avg    # smalls are partial alias fits    # bigs fits and remainder thrown to smalls    smalls = ((i, w/avg) for (i, w) in enumerate(weights) if w < avg)    bigs = ((i, w/avg) for (i, w) in enumerate(weights) if w >= avg)    # [(0, 0.6), (1, 0.9)] --> small weight avgs    # [(2, 1.5)] --> big weight avgs    #cycle through small and big partition generator    small, big = next(smalls, None), next(bigs, None)    # if elems are not exhauted kindly loop    while big and small:        # put the partial elem in aliases        aliases[small[0]] = (small[1], big[0])        print("alias lookup table transformation", aliases)        # big = i of big , weight_b - (1 - weight_s)        # big = (0, (1.5 - (1 - 0.6))        big = (big[0], big[1] - (1-small[1]))        if big[1] < 1:            print("large weight is < 1 skip", aliases)            small = big            big = next(bigs, None)        else:            small = next(smalls, None)    print("alias table generated", aliases)    # SELECTION STAGE    def weighted_random():        r = random.random() * N        i = int(r)        # choose a random probability        # from the alias table        (odds, alias) = aliases[i]        print("what are the odds of", r)        # return the larger probability if        # it's odds are higher else        # return the smaller probability index        return alias if (r-i) > odds else i    return weighted_random()# single trial selectionprint(f"experiment trial", prepare_aliased_randomizer([ 0.2, 0.3, 0.]))# Repeat the experiment to observe bias factorfor e in range(10):  print(f"experiment trial", prepare_aliased_randomizer([0.2, 0.3, 0.5]))

PHEWW!!!!!!! there's alot going on. Let's break it down. There are two major steps, once you understand how the table is generated(that's most of the work) the lookup naturally follows.

This paragraph is interesting but ultimately a digression, you can choose to skip over it!
How exactly does a computer remember who Sam is? Well computers are a little like the human brain, people have complex organic circuitry of neurons that consolidate sensory observation and from that recognizes patterns and form "thoughts", the only difference is we have to tell the computer what pattern it should follow and it doesn't necessary form a thought, simply a result, and we do not need such instructions?(dna? nope, it's a philosophical question of not only intelligence but implicitly will :) does existence precede essence? I like to think so. Shameless plug, @ me on twitter if you're interested in this topic).

Table generation

To create the table, the original weight distribution is sliced into pieces(partitions)! Quite literally! Suppose we know how many weights we need, in this case 3, that's the limit. To find P(0.2, 0.3, 0.5) We create an empty alias, with a data structure, with the exact same length [partition, partition, partition] and a probability [partition, partition, partition]. Each partition will contain exactly P(1/3). Here's the interesting part, we then scale the probabilities to the factor of our limit and separate it into chunks less than 1 and greater than or equal to 1. Why?

Here is the rationale:

1. If the weight fraction is greater than one or equal to it, it will fill a single alias partition(fit exactly if equal else) with some change! which will be sent to an empty partition.

2. If the weight is less than one, then it will be smaller than the allocated partition and it will be expecting a buddy :)

Each partition, will hold a new distribution perfectly and why this is mathematically correct.. is kinda funny tbh. I won't discuss the formal proof, but the idea is that at any point in time, the sum of the elements in the distribution is always proportionate to the original weight.

defmodule Probability do  def bias_probability(weights) do    # Initialization    n = Enum.count(weights)    prepare_alias_table(weights, n)  end  defp prepare_alias_table(weights, n) do    alias_table = Enum.map(1..n, fn _ -> {0, nil} end)    prob = Enum.map(1..n, fn _ -> nil end)    # create work lists    scaled_weight = scale_probability(weights, n)    small =      scaled_weight      |> Enum.filter(fn {_, w} -> w < 1 end)    large =      scaled_weight      |> Enum.filter(fn {_, w} -> w >= 1 end)    # recursively create table (TCO optimized)    transform_alias(small, large, alias_table, prob)  end  # Base case when small and large are empty  defp transform_alias([], [], _, prob), do: prob  defp transform_alias(small = [], [_g = {i, _} | tail], alias_table, prob) do    # Remove the first element from large call it g, Set prob[g] = 1    transform_alias(      small,      tail,      alias_table,      List.replace_at(prob, i, 1)    )  end  defp transform_alias([_l = {i, _} | tail], large = [], alias_table, prob) do    # (clause will trigger due to numerical instability)    # Remove the first element from Small, call it l    transform_alias(      tail,      large,      alias_table,      List.replace_at(prob, i, 1)    )  end  defp transform_alias(         [{index_l, weight_l} | tail_s],         [_g = {index_g, weight_g} | tail_l],         alias_table,         prob       ) do    # Remove the first element from small call it l    # Remove the first element from large call it g    # Pg := (pg + pl) - 1 (numerical stability :) )    new_weight_g = (weight_g + weight_l) - 1    # if Pg < 1 add g to small    if new_weight_g < 1 do      transform_alias(        [{index_g, new_weight_g} | tail_s],        tail_l,        List.replace_at(alias_table, index_l, weight_g),        List.replace_at(prob, index_l, weight_l)      )      # else Pg >= 1 add g to large    else      transform_alias(        tail_s,        [{index_g, new_weight_g} | tail_l],        List.replace_at(alias_table, index_l, weight_g),        List.replace_at(prob, index_l, weight_l)      )    end  end  # HELPER  defp scale_probability(probs, n) do    0..n    |> Enum.zip(probs)    |> Stream.map(fn {i, w} -> {i, w * n} end)  endend

Now that we can generate an alias table, this implementation will also include a local cache using Erlang Term Storage to hold the alias table inside another table! lol throughout the lifecycle of the application once it has begun to run. I'm assuming once a difficulty level is selected random rooms need to be populated but the probability of that room occuring needs to be computed once.

  def bias_probability(weights) do    # Initialization    # does the probability exist in memory?    current_probability = try do      [weights: cached_probs] = :ets.lookup(:weight_probability, :weights)      cached_probs    rescue      ArgumentError -> cache(weights)    end    current_probability    #|> weighted_random # to be implemented next  end  defp cache(weights) do    n = Enum.count(weights)    :ets.new(:weight_probability, [:duplicate_bag, :private, :named_table])    :ets.insert(:weight_probability, {:weights, prepare_alias_table(weights, n)})    [weights: cached_probs] = :ets.lookup(:weight_probability, :weights)    cached_probs  end

Look up

Finally we can generate a random probability and select an index to return from it!

  # GENERATION  defp weighted_random(aliased_table, n) do    # Generate a fair random distro in a range    # from n and call it i.    r = Enum.random(0..1000)/1000 * n    # random choice P(1/3)    i = floor(r) # 0, 1 , 2    prob = aliased_table    {:ok, odd} = Enum.fetch(prob, i)    # partial fit    if (r - i) > (odd) do      # which piece of what goes where      bias = prob      |> Enum.with_index()      |> Stream.filter(fn {p, _} -> p == 1 end)      |> Enum.random()      {_, i} = bias      i    else       i    end  end

Here's the whole implementation in a Github gist

GOTCHAS

There are quite a few waiting for you if you desire to copy/paste this implementation, for example I ignored float point precision(somewhat) instead I could have chosen to seed the values and my probabilities are always weight fractions of 1, or the implicit assumption of never returning 0, or 1 weights in the linear search. Please read the more indepth blogs, they're invaluable. For my use case though, it could probably scale to thousands of rooms effortlessly Thanks for reading!!! I'm always looking for feedback, have a suggestion? found a bug? Know an interesting algorithm? Maybe you just wanna leave kind words, I'm always happy to converse :)