How to simulate a loaded dice in Igor Pro?

Share icon

Hi. I am looking for a way to generate random states with non-equal probabilities. A simple example to illustrate this would be a loaded dice, where  outcome of 1 is 2p/3, the probability of obtaining 2, 3, 4 or 5 is p each, and the probability of obtaining 6 is 3p/2. My approach so far has been to pick random numbers from a wave that would have 6x100 points with the first 66 points equal 1, next 100 equal 2, and so on, with the last 137 points equal to 6. This code shows how a loaded coin would work:

Function return_prob(Probability)
    Variable Probability
    Make/O/N=100/FREE ProbW=1*(x<Probability*100)+0*(x>=Probability*100)
    Return ProbW[50+enoise(50)]
End

However this is computationally very expensive and takes too long to be practical, as I need to generate large sets of data. Also has low resolution, as sometimes the probabilities are at the 7th digit, and making a wave with Nx10^7 only makes things worse.

Is there a function in Igor Pro similar to the one in R that can do this faster/in one shot?  This is how that code would look like:

myprobs <- c(2/3 * p, p, p, p, p, 3/2*p)
x <- sample(1:6, size = 100, replace = TRUE, prob = myprobs)

how about

Function return_prob(Probability)
    Variable Probability
    Return (enoise(0.5)+0.5)<Probability
End

 

March 10, 2021 at 12:39 pm - Permalink

A roll of the die is floor(enoise(3) + 4), so for p(6) > 1/6, try this:

function RollTheDie(p6)
    variable p6
    return floor(min(6, enoise(3) + 3 + 6*p6))
end

 

March 10, 2021 at 11:10 pm - Permalink

In your example, IVP, the probabilities of the loaded dice do not add up to 100%. Is that intentional? If the probability of rolling a 6 is 3/2*1/6 I would expect the probability of rolling a 1 to be 1/2*1/6. This is assuming that the only effect of the loaded dice is to convert a given fraction of 1s into 6s.

To make that clear, I would write the code like this:

Function RollLoadedDice(LoadEfficiency)
//  Rolls a loaded six-sided dice, where a roll of 1 is converted to a roll of 6 with a probability of LoadEfficiency
Variable LoadEfficiency
 
    //  Generates a random value between 0 and 6
    Variable Roll=enoise(3)+3
 
    //  Converts a roll of 1 to a 6 with a probability of LoadEfficiency
    if (Roll<LoadEfficiency)
        Roll=6
        
    //  Rounds up the generated random number to the closest integer, producing rolls of 1, 2, 3, 4, 5 and 6
    else
        Roll=Ceil(Roll)
    endif
 
    //  Returns the rolled value
    Return Roll
end
​

 

March 11, 2021 at 03:04 am - Permalink

or do you need something more general where you have a list of events, each with different probabilities?

event 1      0.10
event 2      0.20
event 3      0.20
event 4      0.40
event 5      0.90
etc...

March 11, 2021 at 03:15 am - Permalink

This is how I would handle the more generalized case. I haven't tested the histogram. You'll want to make sure to do that before you use the code.

Function RandomEventTest()
 
    //  Creates a not-normalised probability wave with 7 events, each with a given probability
    Make/O root:MyProbabilityWave={0.01, 0.05, 0.05, 0.20, 0.30, 0.90, 0.90}
    Wave MyProbabilityWave=root:MyProbabilityWave
 
    //  Creates a list of random events based on the probabilities in MyProbabilityWave
    Make/O/N=100 root:ListOfRandomEvents=RandomEventGenerator(MyProbabilityWave)
    Wave ListOfRandomEvents=root:ListOfRandomEvents
    
    //  Displays the probability wave and the list of random events
    Edit MyProbabilityWave, ListOfRandomEvents
end
 
 
 
 
Function RandomEventGenerator(ProbabilityWave)
//  Selects a random event, based on a list of probabilities for each event
Wave ProbabilityWave
 
    //  Calculates the sum off all probabilities in the wave. For a normalised wave this should be 1. This step is included just in case it is not
    Variable SumOfAllProbabilities=sum(ProbabilityWave)
    
    //  Generates the random number used to select an event
    Variable RandomNumber=enoise(SumOfAllProbabilities/2)+SumOfAllProbabilities/2
    
    //  Integrates the values in the probability wave. This is needed to select an event with the random number
    Duplicate/O/FREE ProbabilityWave IntegratedProbabilityWave
    Integrate ProbabilityWave /D=IntegratedProbabilityWave
 
    //  Counts through the events one at a time and selects the event where the combined probability have passed the random number
    Variable NumberOfEvents=NumPnts(IntegratedProbabilityWave)
    Variable Event=0
    for (Event=0; (RandomNumber>IntegratedProbabilityWave[Event]) && (Event<NumberOfEvents); Event+=1)
    endfor
    
    //  Returns the selected event
    Return Event
end

 

March 11, 2021 at 04:02 am - Permalink