{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using Plots, Statistics, CSV, DataFrames, Printf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Test of a Random Number Generator Using Random Walks\n", "### Gabe Schumm\n", "### Due September 29th, 2021" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Random Number Geneator and Associated Functions" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "P_n (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function rand_int(r)\n", " if typeof(r) != UInt64\n", " r = convert(UInt64, r)\n", " end\n", " \n", " a::UInt64 = 2862933555777941757\n", " c::UInt64 = 1013904243\n", " \n", " return (a*r) + c\n", "end\n", "\n", "function get_bit(r, i::Int64)\n", " if typeof(r) != UInt64\n", " r = convert(UInt64, r)\n", " end\n", " return convert(Int64, (r >>> (i-1)) &1)\n", "end\n", "\n", "\n", "#Stirling's approximation for log(y!)\n", "#Only used when y > 20\n", "function log_fact_appox(y::Int64)\n", " return y *log(y) - y + log(y*(1+(4*y)*(1+2*y)))/6 + log(pi)/2\n", "end\n", "\n", "\n", "function P_n(x::Int64, n::Int64)\n", " \n", " log_n_fact = log_fact_appox(n) #n>20 for this assignment\n", " \n", " np = (n+x) ÷ 2\n", " nm = (n-x) ÷ 2\n", " \n", " if np > 20\n", " log_np_fact = log_fact_appox(np)\n", " else\n", " log_np_fact = log(factorial(np))\n", " end\n", " \n", " if nm > 20\n", " log_nm_fact = log_fact_appox(nm)\n", " else\n", " log_nm_fact = log(factorial(nm))\n", " end\n", " \n", " log_P = log_n_fact - (n * log(2)) - log_np_fact -log_nm_fact\n", " \n", " return exp(log_P)\n", "end\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "random_walk (generic function with 1 method)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function random_walk(N_w::Int64, n::Int64, r_0::UInt64)\n", " #N_w is number of random walk simulations to run\n", " #n is number of steps in each random walk simulation\n", " #r_0 is random UInt64 seed\n", " \n", " counts_array = zeros(64, 2*n + 1) #histogram of final walk distance for each bit\n", " # one row for each bit, 2*n + 1 columns for each\n", " #final distance value (odd and even)\n", " r = r_0\n", " for w =1:N_w\n", " walk = zeros(64) #random walk distance for each of 64 bits\n", "\n", " for i=1:n\n", " r = rand_int(r)\n", " walk += 2 .* get_bit.(r, 1:64) .- 1 #step added for each bit \n", " end\n", "\n", " for b=1:64\n", " x_final = convert(Int64, walk[b] + n + 1) #final distance converted to index of counts_array\n", " counts_array[b, x_final] +=1 #final distance count incremented for each bit \n", " end\n", " end\n", " \n", " delta = zeros(64, 2*n + 1) #array for delta value of each bit and final walk distance\n", " \n", " for b=1:64\n", " delta[b, :] = (counts_array[b, :]/N_w) .- P_n.(-n:n, n) #for each bit, calculate delta\n", " #based on counts array\n", " end\n", " \n", " delta2 = sum(delta[:, 1:2:end] .^2, dims=2) #sum delta^2 along all walk lengths, for each bit\n", " #only include even x values, assuming n is even\n", " return counts_array, sqrt.(delta2) \n", " \n", "end\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "run_random_walk (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function run_random_walk(input)\n", " n_params = size(input)[1] #number of N_w values in input\n", " \n", " N_w_delta = zeros(64, n_params) #sqrt(N_w) * delta_b vs. bit number for each simulation\n", " \n", " #DataFrames containing results of each simulation, for each bit (64 df's)\n", " #x = final position, only even x values, assuming n is even\n", " #P_exact = probability of final position from binomial dist., takes x value and n value\n", " dist_dfs = [DataFrame(x = -100:2:100, P_exact = P_n.(-100:2:100, input[1, \"n\"], )) for i =1:64]\n", " \n", " for i=1:n_params\n", " N_w = input[i, \"N_w\"] #number of random walk simulations to run\n", " \n", " counts_array, delta = random_walk(N_w, input[i, \"n\"], input[i, \"r_0\"]) #run simulation with given paramters\n", " \n", " N_w_delta[:, i] = delta * sqrt(N_w) #calculate sqrt(N_w) * delta_b and store value\n", " \n", " for b=1:64\n", " #column for probability of final position for given number of simulations (N_w)\n", " dist_dfs[b][!, string(\"P_\", N_w)] = counts_array[b,:][1:2:end] / N_w \n", " end\n", " end\n", " \n", " return dist_dfs, N_w_delta\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation Parameters " ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "n_values = 100 * ones(Int64, 6) #n=100 for all 6 runs\n", "N_w_values = 10 .^ (2:7) #Number of random walk simulations run, 10^2 - 10^7\n", "r_0_values = rand(UInt64, 6) #random seed for each run\n", "\n", "input = DataFrame(n = n_values, N_w = N_w_values, r_0 = r_0_values); #organize inputs in DataFrame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run Simulation" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

9 rows × 8 columns

xP_exactP_100P_1000P_10000P_100000P_1000000P_10000000
Int64Float64Float64Float64Float64Float64Float64Float64
1-80.05795840.060.0590.06120.058510.0583970.0578962
2-60.06659050.070.0530.06550.066280.0667630.0666624
3-40.0735270.080.0830.07240.073070.0733710.0735431
4-20.07802870.030.0820.07950.077050.0780730.0781397
500.07958920.080.0790.07620.07920.0792930.0797125
620.07802870.070.0820.07560.078420.0784750.0781047
740.0735270.040.0890.07570.074480.0734330.0735032
860.06659050.030.0650.06720.06740.0664320.0665313
980.05795840.040.050.05880.057250.0579260.0578475
" ], "text/latex": [ "\\begin{tabular}{r|cccccccc}\n", "\t& x & P\\_exact & P\\_100 & P\\_1000 & P\\_10000 & P\\_100000 & P\\_1000000 & P\\_10000000\\\\\n", "\t\\hline\n", "\t& Int64 & Float64 & Float64 & Float64 & Float64 & Float64 & Float64 & Float64\\\\\n", "\t\\hline\n", "\t1 & -8 & 0.0579584 & 0.06 & 0.059 & 0.0612 & 0.05851 & 0.058397 & 0.0578962 \\\\\n", "\t2 & -6 & 0.0665905 & 0.07 & 0.053 & 0.0655 & 0.06628 & 0.066763 & 0.0666624 \\\\\n", "\t3 & -4 & 0.073527 & 0.08 & 0.083 & 0.0724 & 0.07307 & 0.073371 & 0.0735431 \\\\\n", "\t4 & -2 & 0.0780287 & 0.03 & 0.082 & 0.0795 & 0.07705 & 0.078073 & 0.0781397 \\\\\n", "\t5 & 0 & 0.0795892 & 0.08 & 0.079 & 0.0762 & 0.0792 & 0.079293 & 0.0797125 \\\\\n", "\t6 & 2 & 0.0780287 & 0.07 & 0.082 & 0.0756 & 0.07842 & 0.078475 & 0.0781047 \\\\\n", "\t7 & 4 & 0.073527 & 0.04 & 0.089 & 0.0757 & 0.07448 & 0.073433 & 0.0735032 \\\\\n", "\t8 & 6 & 0.0665905 & 0.03 & 0.065 & 0.0672 & 0.0674 & 0.066432 & 0.0665313 \\\\\n", "\t9 & 8 & 0.0579584 & 0.04 & 0.05 & 0.0588 & 0.05725 & 0.057926 & 0.0578475 \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m9×8 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m x \u001b[0m\u001b[1m P_exact \u001b[0m\u001b[1m P_100 \u001b[0m\u001b[1m P_1000 \u001b[0m\u001b[1m P_10000 \u001b[0m\u001b[1m P_100000 \u001b[0m\u001b[1m P_1000000 \u001b[0m\u001b[1m P_100\u001b[0m ⋯\n", "\u001b[1m \u001b[0m│\u001b[90m Int64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float\u001b[0m ⋯\n", "─────┼──────────────────────────────────────────────────────────────────────────\n", " 1 │ -8 0.0579584 0.06 0.059 0.0612 0.05851 0.058397 0.05 ⋯\n", " 2 │ -6 0.0665905 0.07 0.053 0.0655 0.06628 0.066763 0.06\n", " 3 │ -4 0.073527 0.08 0.083 0.0724 0.07307 0.073371 0.07\n", " 4 │ -2 0.0780287 0.03 0.082 0.0795 0.07705 0.078073 0.07\n", " 5 │ 0 0.0795892 0.08 0.079 0.0762 0.0792 0.079293 0.07 ⋯\n", " 6 │ 2 0.0780287 0.07 0.082 0.0756 0.07842 0.078475 0.07\n", " 7 │ 4 0.073527 0.04 0.089 0.0757 0.07448 0.073433 0.07\n", " 8 │ 6 0.0665905 0.03 0.065 0.0672 0.0674 0.066432 0.06\n", " 9 │ 8 0.0579584 0.04 0.05 0.0588 0.05725 0.057926 0.05 ⋯\n", "\u001b[36m 1 column omitted\u001b[0m" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#first(dist_dfs[1], 5)\n", "dist_dfs[64][47:55,:]" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

5 rows × 6 columns

100100010000100000100000010000000
Float64Float64Float64Float64Float64Float64
19.4719129.952894.7191299.528947.1912995.28
29.4719129.952894.7191299.528947.1912995.28
36.3268620.23194.7191299.528632.6862995.28
43.9343720.007363.2686200.073523.1932000.73
54.6383410.977852.3193165.448462.041654.48
" ], "text/latex": [ "\\begin{tabular}{r|cccccc}\n", "\t& 100 & 1000 & 10000 & 100000 & 1000000 & 10000000\\\\\n", "\t\\hline\n", "\t& Float64 & Float64 & Float64 & Float64 & Float64 & Float64\\\\\n", "\t\\hline\n", "\t1 & 9.47191 & 29.9528 & 94.7191 & 299.528 & 947.191 & 2995.28 \\\\\n", "\t2 & 9.47191 & 29.9528 & 94.7191 & 299.528 & 947.191 & 2995.28 \\\\\n", "\t3 & 6.32686 & 20.231 & 94.7191 & 299.528 & 632.686 & 2995.28 \\\\\n", "\t4 & 3.93437 & 20.0073 & 63.2686 & 200.073 & 523.193 & 2000.73 \\\\\n", "\t5 & 4.63834 & 10.9778 & 52.3193 & 165.448 & 462.04 & 1654.48 \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "\u001b[1m5×6 DataFrame\u001b[0m\n", "\u001b[1m Row \u001b[0m│\u001b[1m 100 \u001b[0m\u001b[1m 1000 \u001b[0m\u001b[1m 10000 \u001b[0m\u001b[1m 100000 \u001b[0m\u001b[1m 1000000 \u001b[0m\u001b[1m 10000000 \u001b[0m\n", "\u001b[1m \u001b[0m│\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\n", "─────┼───────────────────────────────────────────────────────\n", " 1 │ 9.47191 29.9528 94.7191 299.528 947.191 2995.28\n", " 2 │ 9.47191 29.9528 94.7191 299.528 947.191 2995.28\n", " 3 │ 6.32686 20.231 94.7191 299.528 632.686 2995.28\n", " 4 │ 3.93437 20.0073 63.2686 200.073 523.193 2000.73\n", " 5 │ 4.63834 10.9778 52.3193 165.448 462.04 1654.48" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = DataFrame(N_w_delta, :auto)\n", "d = rename(df, Dict([names(df)[i] => string(N_w_values[i]) for i=1:6]))\n", "first(d, 5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Export Simulation results" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "CSV.write(\"d.csv\", d); " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "for b=0:63\n", " filename = @sprintf(\"p%02i.csv\",b)\n", " CSV.write(string(\"bit_prob_dists/\", filename), dist_dfs[b+1])\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exploring Periodicity of Each Bit" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [], "source": [ "r = 502\n", "sequence = []\n", "\n", "b = 5 #bit of interest\n", "p = 2^b #periodicity of bit of interest\n", "\n", "for i=1:100\n", " r = rand_int(r)\n", " push!(sequence, get_bit(r, b))\n", "end" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "true\n", "0\n" ] } ], "source": [ "p=2^b\n", "println(sequence[1:p] == sequence[1+p:2*p]) #first sequence of p bits = next sequence of p bits \n", "\n", "println(sum(2 .* sequence[1:p] .- 1)) #sum of cycle of steps = 0" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 & 110 & 2 & 6\n", "1 & 001 & 1 & 1\n", "2 & 000 & 0 & 0\n", "3 & 011 & 3 & 3\n", "4 & 010 & 2 & 2\n", "5 & 101 & 1 & 5\n", "6 & 100 & 0 & 4\n", "7 & 111 & 3 & 7\n", "8 & 110 & 2 & 6\n", "9 & 001 & 1 & 1\n", "10 & 000 & 0 & 0\n", "11 & 011 & 3 & 3\n", "12 & 010 & 2 & 2\n", "13 & 101 & 1 & 5\n", "14 & 100 & 0 & 4\n" ] } ], "source": [ "r=502\n", "for i =1:15\n", " #print iteration number, first three bits of string, r mod 4, and r mod 8\n", " println(i-1, \" & \", string(bitstring(r)[62:end]), \" & \", r % 4, \" & \", r % 8)\n", " r = rand_int(r)\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.2", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }