# How to simulate a loaded dice in Igor Pro? 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)

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

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