Let’s start with a *very* brief Mathematica primer.

(* Mathematica comments appear between parentheses and asterisks like this. *)

(* Here is an example of a function in Mathematica. *)

(* f1 is a function of a single input (x). *)
f1[x_]:= Module[
{localvariable1, localvariable2},
localvariable1 = 2 * x;
localvariable2 = localvariable1 + 5;
(* If the line does not end in a semicolon, then it will be the output of the function *)
localvariable2
];

(* See how this function behaves on an input. *)
f1[5]

(* You can loop over a list of possibilities with the Mathematica function "Table" *)

(* I read this as: "data is a list of i and f1[i], with i ranging from 1 to 11 in increments of 0.4" *)

(* If there were no semicolon on the end, Mathematica would print out the result for you to see. *)
data = Table[
{i, f1[i]},
{i, 1, 11, 0.4}
];

(* You can plot that data using the function ListPlot *)
ListPlot[data]

(* I like to make my plots a little fancier *)
ListPlot[
data,
Frame->True,
AspectRatio->1,
ImageSize->200,
PlotMarkers->{Automatic, 8},
FrameLabel->{"x coordinate", "y coordinate"}
]

Now let’s deal with the 1D random walker (p = q = 0.5)

T = 100; (* Number of steps or time *)
p = 0.5; (* Probability of stepping to the right *)
q = 1-p;

(* The Mathematica function Table[] creates a list of +1 and -1 for each step, indicating whether we would take a step to the right or to the left.*)
Steps = Table[If[RandomReal[]<p, 1, -1], {step, 1, T}]

(* Adding up all the steps gives the value of X after T steps. (The Mathematica function Total[] sums up all the elements of a list.) *)
Total[Steps]

(* We can repeat this procedure NumberOfWalkers times to generate statistics for a histogram *)
NumberOfWalkers = 10000;
XValuesFromLotsOfWalkers=Table[
Total[Table[If[RandomReal[] < p, 1, -1], {step, 1, T}]],
{WalkerNumber, 1, NumberOfWalkers} 
];

(* Now we use Mathematica's Histogram[] function to plot this distribution, using bins that are 1 unit wide. *)
Histogram[XValuesFromLotsOfWalkers, {1} (* This 1 sets how wide the bins are *)]

(* If I want to control my axis labels as before, I can using the Show[] function. This also lets me combine multiple plots... *)
Show[
(* Here is the plot with the histogram for the sampled walkers *)
Histogram[
XValuesFromLotsOfWalkers, {1}
],
(* Here is a plot with the exact probability distribution, computed as our (Subscript[P, T])(X) *)
ListPlot[
ExactProbabilityDistribution = Table[
{X,NumberOfWalkers*(Factorial[T]/(Factorial[(T+X)/2] * Factorial[(T-X)/2]))(p^((T+X)/2))(q^((T-X)/2))}, {X, -40, 40}
]
],
Frame->True,
AspectRatio->1,
ImageSize->200,
PlotMarkers->{Automatic, 8},
FrameLabel->{"X", 
ToString[Subscript["P", ToString[T]], StandardForm]<>"(X)",
 "1D random walker, p=q"}
]

(* Similarly I can calculate the mean and variance with Mathematica's Mean[] and Variance[] functions. The function N[] is being used to get a numerical approximation (a decimal) rather than having a fractional answer. *)
mean = N[Mean[XValuesFromLotsOfWalkers]]
variance = N[Variance[XValuesFromLotsOfWalkers]]

(* You'll notice that the mean is not exactly zero. This is because we have "statistical error". We have only looked at NumberOfWalkers random walkers. If we increase the number of walkers, we will see that error decrease, but it won't go to zero unless we have an infinite number of random walkers. In the next section we address how to report the statistical error with error bars. *)

One d random walker (p≠q; r=0)

T = 100;
p = 0.65;
q = 1 - p;

NumberOfWalkers = 10000;
XValuesFromLotsOfWalkers=Table[
Total[Table[If[RandomReal[] < p, 1, -1], {step, 1, T}]],
{WalkerNumber, 1, NumberOfWalkers}
];

(* Notice that the mean is now skewed to the right because p > q *)
Show[
Histogram[
XValuesFromLotsOfWalkers, 50
],
Frame->True,
AspectRatio->1,
ImageSize->200,
PlotMarkers->{Automatic, 8},
FrameLabel->{"X",
ToString[Subscript["P", ToString[T]], StandardForm]<>"(X)", "1D random walker, p="<>ToString[p]<>", q="<>ToString[q]}
]
mean = N[Mean[XValuesFromLotsOfWalkers]]

(* To estimate the size of the error bars, we should repeat the procedure at least ten times and report the variance of those repeats. *)

NumberOfWalkersPerTrial = 1000;
NumberOfTrials = 10;

ListOfMeanValuesFromEachTrial=
Table[
XValuesFromLotsOfWalkers=Table[
Total[Table[If[RandomReal[] < p, 1, -1], {step, 1, T}]],
{WalkerNumber, 1, NumberOfWalkersPerTrial}
];
N[Mean[XValuesFromLotsOfWalkers]],
{i, 1, NumberOfTrials}
]

(* The value of <X(Subscript[>, T]) is then ... *)
Mean[ListOfMeanValuesFromEachTrial]
(* plus or minus *)
2StandardDeviation[ListOfMeanValuesFromEachTrial]
(* By making error bars that extend out to 2 times the standard deviation, we expect our error bars to include the right answer (which we computed exactly to be (p-q)T) 95% of the time *)

(* We could repeat this for various choices of T, computing the error bars each time... *)

NumberOfWalkersPerTrial = 1000;
NumberOfTrials = 10;

DataForManyTimes = Table[
(* This table loops choices of T *)

(* For each T, we run NumberOfTrials different trials, each trial consisting of NumberOfWalkersPerTrial random walkers *)

ListOfMeanValuesFromEachTrial=
Table[
(* This loops over NumberOfTrials different trials *)
XValuesFromLotsOfWalkers=Table[
(* This loops over NumberOfWalkersPerTrial random walkers *)
Total[Table[If[RandomReal[] < p, 1, -1], {step, 1, T}]],
{WalkerNumber, 1, NumberOfWalkersPerTrial}
];
N[Mean[XValuesFromLotsOfWalkers]],
{i, 1, NumberOfTrials}
];

(* The loop over T returns the value of T, the average <X(Subscript[>, T]) estimated from sampling the walkers, and the error bars for that estimate *)
{T, Mean[ListOfMeanValuesFromEachTrial], 2StandardDeviation[ListOfMeanValuesFromEachTrial]},

{T, {1, 10, 25, 50, 100}}
]

Needs["ErrorBarPlots`"]
Show[
ErrorListPlot[DataForManyTimes,
PlotMarkers->{Automatic, 12}
],
Plot[(p-q)T, {T, 0, 100},
PlotStyle->{Black, Dashed}
],
Frame->True,
AspectRatio->1,
ImageSize->300,
PlotMarkers->{Automatic, 8},
FrameLabel->{Style["T", FontSize->18],
Style["<X(Subscript[>, T])", FontSize->18], Style["1D random walker, p="<>ToString[p]<>", q="<>ToString[q], FontSize->18]}
]

(* Instead of making a bar plot histogram, I can get the height of the bars with the function HistogramList[] *)
NumberOfWalkers = 10000;
T = 10;
XValuesFromLotsOfWalkers=Table[
Total[Table[If[RandomReal[] < p, 1, -1], {step, 1, T}]],
{WalkerNumber, 1, NumberOfWalkers} 
];

{Xvalues, Xcounts} = HistogramList[XValuesFromLotsOfWalkers, {1}(*Bins of width 1*)]

(* Notice that every other bar has a height of zero, so let's exclude those. If I have a list, I can get elements from i to j counting by k if I do list[[i;;j;;k]]. In other words list[[1;;100;;3]] would give back the first, fourth, seventh, ... elements of the list. Also notice that there is a shorthand for the last item of a list. Since you might not know how long the list is, you sometimes want to just indicate the last item with -1. Or the second to last item with -2.
*)

(* So the following will give a list of the X coordinate and the fraction of times that X coordinate was measured *)
HistogramValuesT10 = Transpose[
{Xvalues[[1;;-2;;2]], Xcounts[[;;;;2]] / NumberOfWalkers}
]

(* You can repeat this with a bigger value of T *)
NumberOfWalkers = 10000;
T = 25;
XValuesFromLotsOfWalkers=Table[
Total[Table[If[RandomReal[] < p, 1, -1], {step, 1, T}]],
{WalkerNumber, 1, NumberOfWalkers} 
];
{Xvalues, Xcounts} = HistogramList[XValuesFromLotsOfWalkers, {1}(*Bins of width 1*)];
HistogramValuesT25 = Transpose[
{Xvalues[[1;;-2;;2]], Xcounts[[;;;;2]] / NumberOfWalkers}
];

(* And plot both together with *)
ListPlot[
(* This next line is the list of things you're plotting. (If you were wanting to do a bunch of choices of T it might be smarter to loop over the choices by using a Table) *)
{HistogramValuesT10, HistogramValuesT25},

(* Below here is setting the style of the plot *)
Joined->True,
PlotStyle->Dashed,
Frame->True,
AspectRatio->1,
ImageSize->300,
Axes->False,
PlotMarkers->{Automatic, 8},
FrameLabel->{Style["X", FontSize->18],
Style["(Subscript[P, T])(X)", FontSize->18], Style["1D random walker, p="<>ToString[p]<>", q="<>ToString[q], FontSize->18]},
PlotLegends->{"T = 10", "T=25"}
]

A note about the p, q, r case

(* The earlier simulations were making a list of all the steps and marking them down as +1 or -1. Now we need to include the possibility of 0. Remember that this is what we were using when there were only two options: *)

(* The Mathematica function Table[] creates a list of +1 and -1 for each step, indicating whether we would take a step to the right or to the left.*)
Steps = Table[If[RandomReal[]<p, 1, -1], {step, 1, T}]

(* Now if I were trying to pick from three options (p, q, or r), I could do it like: *)
RandomNumberIDrew = RandomReal[];
value = If[RandomNumberIDrew < p, 1, If[RandomNumberIDrew < (p + q), -1, 0] ]
(* The random number is drawn between 0 and 1. This is first checking if the random number is less than p (which occurs with probability p). If so, value is being set equal to 1. If not, the next If statement is used to check if the random number is less than (p + q). We already knew that the number was bigger than p (otherwise we wouldn't have made it to this step, so it will be less than (p + q) with probability q. If that occurs, value is set equal to -1 (a step to the left). Finally, if neither happens (probability r = 1 - p - q), we take no step by setting value equal to zero.