Distribution Function Sampling

http://arxiv.org/PS_cache/cond-mat/pdf/0104/0104215.pdf

Load needed packages only once

<< Statistics`DataManipulation`
<< Graphics`Graphics`

Uncomment and use this in place of data statements below if data is in a file

Central Limit Theorem sampling of the Gaussian Distribution Function

For Y = the sum of n random numbers in the range (0-1), Y tends to a Gaussian distribution having mean n/2 and variance n/12. Thus q = ((Y - (n/2))/(Sqrt[n/12]))tends to a Gausian distribution having mean 0 and variance 1.

Print["n = ", n  =  12]
data = 
        Table[(((∑_(j = 1)^n) Random[] - (n/2))/(Sqrt[(n/12)])), {i, 1, 
            10000}];
ListPlot[data];

Bin Plots

Print["Minimum x value = ", minX = Min[data]];
Print["Maximum x value = ", 
        maxX = Max[data]];
Print["Bin size = ",  
        binSize = N[(maxX - minX)/50]];
Print["Bin counts = ", 
        binCount = 
          BinCounts[
            data, {minX, maxX, binSize}]];
Histogram[binCount, FrequencyData -> True];
points = 
        Table[{Extract[Range[minX, maxX, binSize], i], 
            Extract[binCount, i]}, {i, 1, Length[binCount]}];
freqPlot = 
        ListPlot[points, PlotStyle -> PointSize[0.04], 
          PlotRange -> {{0, 1}, {0, Max[binCount]}}, 
          DisplayFunction -> Identity];
pointPlot = 
        ListPlot[points, PlotStyle -> PointSize[0.9], 
          PlotJoined -> True, 
          PlotRange -> {{minX, maxX}, {0, Max[binCount]}}, 
          DisplayFunction -> Identity];

    
\(\(\(DisplayTogetherArray[Show[pointPlot, freqPlot], 
        DisplayFunction -> $DisplayFunction];\)\(
\)
    \)\)

Metropolis Monte Carlo sampling of an exponential distribution f[x]

Clear[x, f]
α = 1
f[x_] := Exp[-α*x]
Print["Distribution function: ", f[x]]
Plot[Evaluate[f[x]], {x, 0, 5}, AspectRatio -> 1, 
        Frame -> True];

Make a Markov chain of length n over nstates states

nstates = 10;
n  =  
        10000;
data  =  {0};
x = 
        Random[Integer, {1, 
            nstates}];
Print["Reasonable first \
choice is at max of f => x = ", x = f[0]];
Do[(
xt = 
            Random[Integer, {1, nstates}]; 
w = 
            f[xt]/f[x]; 
If[w >= 1,  {x = xt}, 
            If[Random[] <= w, {x = xt}]]; 
data  =  
            Append[data,  x]), {n - 1}]

Bin Plots

Print["Minimum x value = ", minX = Min[data]];
Print["Maximum x value = ", 
        maxX = Max[data]];
Print["Bin size = ",  
        binSize = N[(maxX - minX)/nstates]];
Print["Bin counts = ", 
        binCount = BinCounts[data, {minX, maxX, binSize}]];
Histogram[binCount, 
        FrequencyData -> True];
points = 
        Table[{Extract[Range[minX, maxX, binSize], i], 
            Extract[binCount, i]}, {i, 1, Length[binCount]}];
freqPlot = 
        ListPlot[points, PlotStyle -> PointSize[0.04], 
          PlotRange -> {{0, 1}, {0, Max[binCount]}}, 
          DisplayFunction -> Identity];
pointPlot = 
        ListPlot[points, PlotStyle -> PointSize[0.9], 
          PlotJoined -> True, 
          PlotRange -> {{minX, maxX}, {0, Max[binCount]}}, 
          DisplayFunction -> Identity];

    
\(\(\(DisplayTogetherArray[Show[pointPlot, freqPlot], 
        DisplayFunction -> $DisplayFunction];\)\(
\)
    \)\)

Cumulative Distribution Function F[x] for an exponential distribution

Clear[x, f]
f[x_] := 
      Exp[-α*x] /. α -> 1
F[x_] := (∫_0^x)(f[y] /. α -> 
                1) \[DifferentialD]y
f[x]
Plot[F[x], {x, 0, 5}, AspectRatio -> 1, 
        Frame -> True];

Inverse of Cumulative Distribution Function

data = 
        Table[-Log[1 - Random[]], {i, 1, 
            1000}];
ListPlot[data];

Bin Plots

Print["Minimum x value = ", minX = Min[data]];
Print["Maximum x value = ", 
        maxX = Max[data]];
Print["Bin size = ",  
        binSize = N[(maxX - minX)/20]];
Print["Bin counts = ", 
        binCount = BinCounts[data, {minX, maxX, binSize}]];
Histogram[binCount, 
        FrequencyData -> True];
points = 
        Table[{Extract[Range[minX, maxX, binSize], i], 
            Extract[binCount, i]}, {i, 1, Length[binCount]}];
freqPlot = 
        ListPlot[points, PlotStyle -> PointSize[0.04], 
          PlotRange -> {{0, 1}, {0, Max[binCount]}}, 
          DisplayFunction -> Identity];
pointPlot = 
        ListPlot[points, PlotStyle -> PointSize[0.9], 
          PlotJoined -> True, 
          PlotRange -> {{minX, maxX}, {0, Max[binCount]}}, 
          DisplayFunction -> Identity];

    
\(\(\(DisplayTogetherArray[Show[pointPlot, freqPlot], 
        DisplayFunction -> $DisplayFunction];\)\(
\)
    \)\)

Miscelaneous

Simple 2D random walk

points  =  {{0,  0}};
n  =  
        100;
x = 
        0;
y = 0;
\(Do[
If\ [Random[Integer]\  == \ 0, \ \(x++\), \ \(x--\)]; 

        If\ [Random[Integer]\  == \ 0, \ \(y++\), \ \(y--\)]; 

        points\  = \ Append[points, \ {x, \ y}], 
{n}]\)
ptPlot  = 
        ListPlot[points, 
          PlotStyle -> {PointSize[0.03], RGBColor[0, 0, 1]}, 
          DisplayFunction -> 
            Identity];
stepPlot  = 
        ListPlot[points, PlotJoined -> True, 
          DisplayFunction -> 
            Identity];
Show[{ptPlot, stepPlot}, 
        AspectRatio -> Automatic, 
        DisplayFunction -> $DisplayFunction];

Boneyard

Make a list of random numbers between 1 and 100

n = 
        10;
lst = {0};
Do[lst = Append[lst,  Random[Integer, {1, 100}]], {n}]
lst

Failed attempt to make a list of Monte Carlo selected distribution

nextx[
        x_] := (
xt = 
          Random[Integer, {1, 100}]; 
w = 
          f[xt]/f[x]; 
  
If[w >= 1, Return[x], 
          If[Random[] <= w, Return[xt], 
            Return[x]]]
)

x = Random[Integer, {1, 100}];
data = Table[nextx[x], {i, 1, 10}]

      
data = Delete[data, Table[{i}, {i, 1, n/10}]];
ListPlot[data];

Print["Minimum x value = ", minX = Min[data]];
Print["Maximum x value = ", 
        maxX = Max[data]];
Print["Bin size = ",  
        binSize = N[(maxX - minX)/20]];
binCount = 
        BinCounts[data, {minX, maxX, binSize}];
barPlot = BarChart[binCount];