<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 I found this data set on Kaggle. It is 1 million single digit numbers from 1 to 9, and I was curious if their distribution is uniform. Or said another way, how reasonable is it suggest the numbers in the list were generated uniformly at random? I've recently learned the basics of using the PyMC3's Markov Chain Monte Carlo methods to answer these types of questions with a Baysian framework, and I think the ideas and methods can be used for this. Lets start with the data.
One Million Pseudo-Random Digits
Step1:  First thing to notice is that the digits are "arranged into ten columns for your convenience," as the data set description states. I think it wil be easier to work with as a flattened set of digits.
Step2:  Now, do the individual digit counts make sense supposing these 1 million digitis were drawn uniformly at random?
Step3:  Zero was left out of the million single digits for some reason. The other digits appear with very similar counts varied around
$$\frac{1000000}{9} = \text{~}111111$$
as you'd expect from a list proporting to be uniformly generated. More, the differences are tiny compared to $1000000/9$, so there is obvious problem with saying these digits were drawn from a uniform distribution. But with a Bayesian framework, and PyMC3, I can ask a subtler question
Step5:  Aside
Most of the reasons for my failures involve exceptions. In IPython however, the error output is too large to be a reasonable part of a narrative, at least I find it so. I still want to report back the exceptions I encounter, so heres a quick function to avoid some boiler plate.
Step6:  Fail 1
The first thing I tried was based on the fact that PyMC3 has a Discrete Uniform distribution in it. In the PyMC3 approach, you link data to disributions with the observed keyword, and run an MCMC algorithm to find posterior distribution(s) to fit your data. The million digits ostensibly is given by a discrete, uniform distribution over the digits 1 to 9, so the PyMC3 model should link the given data with that distribution. (the CD-P book above explains this code)
Step7:  Seems small for being the whole model, but lets run it using the Metroplois algorithm.
Step8:  Fail! The error message gives a hint at what went wrong though. The one variable in the model, digits, is fixed to the data, so there is no posterior distribution being calculated. The "one array to concatenate" is the one array representing a variable we're calculating a posterior distribution for. There are none, so... fail!
Fail 2 & 3
To fix fail 1, I came up with the idea of taking another uniform distribution over 1 to 9, and look at the differences to the data. Not sure why I thought this might be a good idea; I was spit balling. But it failed in another way.
Step9:  It seems likely that I got this memory error because the million data points is too large, so I'll try a subset of the .
Step10:  It runs this time, but the output is odd, sometimes. On occasion, it outputs progress bar over several lines, but not always. Regardless lets see what the "diff" posterior looks like.
Step11:  I take a burned trace to remove the typically wild fluctuations at the beginning of the MC algorithm where the priors are the least constrained, and then I try to chart the results using PyMC traceplot.
Step12:  The difference graph doesn't complete, even after many hours. Something is wrong with this formulation, and I'm going to try a different approach. 
SUCCESS
The idea this time is to build a model where I estimate the probabilities of the digits themselves, and use the PyMC3 categorical distribution with those probabilities as the prior data generation process. Then I can compare the entropy of the estimated distribution fitting the digit data to the entropy of a 9 element uniform distribution, $U(9)$, which is
$$H(U(9)) = \sum_i p_i \log p_i^{-1} = \sum_i \frac{\log 9}{9} = \log 9.$$
I should mention that this is the maximum entropy of a distribution over 9 elements, as this comes up below.
Noe, I could just calculate the entropy of the given set of digits. The entropy estimate of finite set is
$$H(S) = -\sum_i \frac{n_i}{N}\text{log}\left( \frac{n_i}{N} \right)$$
where $n_i$ is the count of the $i$th unique element and $N$ is the size of $S$. The estimate for the of digits is...
Step13:  which is the same to 5 decimales places. Now that I am considering entropy, that proximity is enough to say that in terms of frequency the 1 million digits are generted approximately uniformly, but lets throw PyMC3 at the problem anyways to practice formulating the problem.
To start, I need 9 prior probabilities, one for the probability of drawing each of the nine digits, to define other 9 element distributions.  Each prior is a probability valued real number, so a good prior is a uniform distribution from 0 to 1. I'll use a test value of $1/9$ to initialize $p_1$, so the algorithm will converges more quickly assuming I am correct that $p_1$ should be close to $1/9$.
Step14:  There is one constraint, that the 9 probabilities sum to 1, so $p_2$ can't have the same prior since that would allow $p_1 + p_2$ to be greater than $1$ in some cases. Given $p_1$, $p_2$ can range from $0$ to $1-p_1$. A PyMC prior can be defined so it is dependent on other variables, like $p_1$.
Step15:  The next 6 are similar, each depending on the previously defined probabilities.
Step16:  Finally, $p_9$ is deterministic given the other 8, as the sum must be 1.
Step17:  In the failed models above, I implicitly assumed the probabilities were 1/9 by using the discrete uniform distribution. Here, I let the probabilities wander to match the categorical data. There is one small hitch with using the PyMC categorical distribution; it indexes from 0 to $n$. But the digits in the digit data are from 1 to 9. To get around that, I added $p_0$ with a 0 probability so that the digit 0 has a zero probability instead of something else like subtracting 1 from each data point.
Step18:  Next, I want to calculate the entropy, so that we can see if it comes out near $\log 9$. I'm first going to create another tensor of the 9 probabilities, and calculate the entropy from elementwise operations.
Step19:  But instead of charting the entropy itself, I'll chart $e^H$, which we expect to be close to $e^{\log 9} = 9$, so that the values coming out are close to a whole number. I find this makes the charts easier to read, and I don't think it adds much overhead to the algorithm.
Step20:  Finally, I run the MC algorithm over on the model using the algorithm and sample size I've observed in the learning materials.
<ASSISTANT_TASK:>
Python Code:
import pandas as pd
the_digits = pd.read_csv('data.csv')
the_digits.head()
Explanation: I found this data set on Kaggle. It is 1 million single digit numbers from 1 to 9, and I was curious if their distribution is uniform. Or said another way, how reasonable is it suggest the numbers in the list were generated uniformly at random? I've recently learned the basics of using the PyMC3's Markov Chain Monte Carlo methods to answer these types of questions with a Baysian framework, and I think the ideas and methods can be used for this. Lets start with the data.
One Million Pseudo-Random Digits
End of explanation
the_digits = the_digits.values.flatten()
Explanation: First thing to notice is that the digits are "arranged into ten columns for your convenience," as the data set description states. I think it wil be easier to work with as a flattened set of digits.
End of explanation
from collections import Counter
digit_counts = Counter(the_digits)
digit_counts
for digit in digit_counts:
    print(round(digit_counts[digit] - 1000000/9))
Explanation: Now, do the individual digit counts make sense supposing these 1 million digitis were drawn uniformly at random?
End of explanation
import pymc3 as pm
import theano.tensor as tt
Explanation: Zero was left out of the million single digits for some reason. The other digits appear with very similar counts varied around
$$\frac{1000000}{9} = \text{~}111111$$
as you'd expect from a list proporting to be uniformly generated. More, the differences are tiny compared to $1000000/9$, so there is obvious problem with saying these digits were drawn from a uniform distribution. But with a Bayesian framework, and PyMC3, I can ask a subtler question: how likely is it that this sample was drawn from a uniform distribution over the 9 digits? There are other distributions over the 9 digits which could have produced the observed data set, and one of those might be a better candidate for generatng distribution of these 1 million digits. 
I'm going attempt to frame thise question in a Bayesian framework using MCMC in PyMC3 to compute the other candidates. 
MCMC in PyMC3
I've recently worked through Cam Davison-Pilon's great book "Probabilistic Programming and Bayesian Methods for Hackers." Well actualy, I worked through the PyMC3 port of that same book by Max Margenot and Thomas Wiecki. I'm going to leave in some of my failed attempts, as seeing the error messages and wrong thinking will be helpful for troubleshooting my problem on other projects. Like always... more imports!
End of explanation
def print_err(err):
    Print the error name and message 
    print(type(err).__name__,':', err)
Explanation: Aside
Most of the reasons for my failures involve exceptions. In IPython however, the error output is too large to be a reasonable part of a narrative, at least I find it so. I still want to report back the exceptions I encounter, so heres a quick function to avoid some boiler plate.
End of explanation
with pm.Model() as fail1:
    digits = pm.DiscreteUniform('bias', 1,9, observed=the_digits)
Explanation: Fail 1
The first thing I tried was based on the fact that PyMC3 has a Discrete Uniform distribution in it. In the PyMC3 approach, you link data to disributions with the observed keyword, and run an MCMC algorithm to find posterior distribution(s) to fit your data. The million digits ostensibly is given by a discrete, uniform distribution over the digits 1 to 9, so the PyMC3 model should link the given data with that distribution. (the CD-P book above explains this code)
End of explanation
with fail1:
    
    try:
        step = pm.Metropolis()
        
    except ValueError as err:
        print_err(err)
    else:
        trace = pm.sample(10000, step=step)
Explanation: Seems small for being the whole model, but lets run it using the Metroplois algorithm.
End of explanation
with pm.Model() as fail2:
    digits = pm.DiscreteUniform('digits', 1,9, observed=the_digits)  # from fail 1
    base = pm.DiscreteUniform('base', 1,9) # the other uniform 
    
    diff = pm.Deterministic('diff', base-digits) # look at the difference
    
    try:
        trace = pm.sample(10000, step=pm.Metropolis())
    except MemoryError as err:
        print_err(err)
Explanation: Fail! The error message gives a hint at what went wrong though. The one variable in the model, digits, is fixed to the data, so there is no posterior distribution being calculated. The "one array to concatenate" is the one array representing a variable we're calculating a posterior distribution for. There are none, so... fail!
Fail 2 & 3
To fix fail 1, I came up with the idea of taking another uniform distribution over 1 to 9, and look at the differences to the data. Not sure why I thought this might be a good idea; I was spit balling. But it failed in another way.
End of explanation
with pm.Model() as fail3:
    base = pm.DiscreteUniform('base', 1,9)
    digits = pm.DiscreteUniform('digits', 1,9, observed=the_digits[:10000]) # use first 10000 points
    diff = pm.Deterministic('diff', base-digits)
    
    step = pm.Metropolis()
    trace = pm.sample(10000, step=step)
Explanation: It seems likely that I got this memory error because the million data points is too large, so I'll try a subset of the .
End of explanation
%matplotlib inline
Explanation: It runs this time, but the output is odd, sometimes. On occasion, it outputs progress bar over several lines, but not always. Regardless lets see what the "diff" posterior looks like.
End of explanation
burned_trace = trace[1000:]  # remove the first 1000 instances, as this is the "burn in period" of the algorithm
pm.plots.traceplot(trace=burned_trace, varnames=['diff']);
Explanation: I take a burned trace to remove the typically wild fluctuations at the beginning of the MC algorithm where the priors are the least constrained, and then I try to chart the results using PyMC traceplot.
End of explanation
from math import log, exp
N = 1000000
entropy = log(9)
ent_est = -sum((n/N) * log(n/N) for n in digit_counts.values())
entropy, ent_est, entropy - ent_est
Explanation: The difference graph doesn't complete, even after many hours. Something is wrong with this formulation, and I'm going to try a different approach. 
SUCCESS
The idea this time is to build a model where I estimate the probabilities of the digits themselves, and use the PyMC3 categorical distribution with those probabilities as the prior data generation process. Then I can compare the entropy of the estimated distribution fitting the digit data to the entropy of a 9 element uniform distribution, $U(9)$, which is
$$H(U(9)) = \sum_i p_i \log p_i^{-1} = \sum_i \frac{\log 9}{9} = \log 9.$$
I should mention that this is the maximum entropy of a distribution over 9 elements, as this comes up below.
Noe, I could just calculate the entropy of the given set of digits. The entropy estimate of finite set is
$$H(S) = -\sum_i \frac{n_i}{N}\text{log}\left( \frac{n_i}{N} \right)$$
where $n_i$ is the count of the $i$th unique element and $N$ is the size of $S$. The estimate for the of digits is...
End of explanation
with pm.Model() as model:
    p1 = pm.Uniform('p1', 0, 1, testval=1/9)
    
Explanation: which is the same to 5 decimales places. Now that I am considering entropy, that proximity is enough to say that in terms of frequency the 1 million digits are generted approximately uniformly, but lets throw PyMC3 at the problem anyways to practice formulating the problem.
To start, I need 9 prior probabilities, one for the probability of drawing each of the nine digits, to define other 9 element distributions.  Each prior is a probability valued real number, so a good prior is a uniform distribution from 0 to 1. I'll use a test value of $1/9$ to initialize $p_1$, so the algorithm will converges more quickly assuming I am correct that $p_1$ should be close to $1/9$.
End of explanation
with model:
    p2 = pm.Uniform('p2', 0, 1-p1, testval=1/9)
    
Explanation: There is one constraint, that the 9 probabilities sum to 1, so $p_2$ can't have the same prior since that would allow $p_1 + p_2$ to be greater than $1$ in some cases. Given $p_1$, $p_2$ can range from $0$ to $1-p_1$. A PyMC prior can be defined so it is dependent on other variables, like $p_1$.
End of explanation
with model:
    p3 = pm.Uniform('p3', 0, 1-p1-p2, testval=1/9)
    p4 = pm.Uniform('p4', 0, 1-p1-p2-p3, testval=1/9)
    p5 = pm.Uniform('p5', 0, 1-p1-p2-p3-p4, testval=1/9)
    p6 = pm.Uniform('p6', 0, 1-p1-p2-p3-p4-p5, testval=1/9)
    p7 = pm.Uniform('p7', 0, 1-p1-p2-p3-p4-p5-p6, testval=1/9)
    p8 = pm.Uniform('p8', 0, 1-p1-p2-p3-p4-p5-p6-p7, testval=1/9)
    
Explanation: The next 6 are similar, each depending on the previously defined probabilities.
End of explanation
with model:
    p9 = pm.Deterministic('p9', 1-p1-p2-p3-p4-p5-p6-p7-p8)
 
Explanation: Finally, $p_9$ is deterministic given the other 8, as the sum must be 1.
End of explanation
with model:
    p_ext = tt.stack([0,p1,p2,p3,p4,p5,p6,p7,p8,p9])
    
    digits = pm.Categorical('digits', p_ext, observed=the_digits)
Explanation: In the failed models above, I implicitly assumed the probabilities were 1/9 by using the discrete uniform distribution. Here, I let the probabilities wander to match the categorical data. There is one small hitch with using the PyMC categorical distribution; it indexes from 0 to $n$. But the digits in the digit data are from 1 to 9. To get around that, I added $p_0$ with a 0 probability so that the digit 0 has a zero probability instead of something else like subtracting 1 from each data point.
End of explanation
with model:    
    p = tt.stack([p1,p2,p3,p4,p5,p6,p7,p8,p9])
    logp = tt.log(p)
    entropy = -tt.sum(tt.mul(logp,p))
Explanation: Next, I want to calculate the entropy, so that we can see if it comes out near $\log 9$. I'm first going to create another tensor of the 9 probabilities, and calculate the entropy from elementwise operations.
End of explanation
with model: 
    e_entropy = pm.Deterministic('e^entropy', tt.exp(entropy))
Explanation: But instead of charting the entropy itself, I'll chart $e^H$, which we expect to be close to $e^{\log 9} = 9$, so that the values coming out are close to a whole number. I find this makes the charts easier to read, and I don't think it adds much overhead to the algorithm.
End of explanation
with model:
    step = pm.Metropolis()
    trace = pm.sample(20000, step=step)
burned_trace = trace[1000:]
pm.plots.traceplot(trace=burned_trace, varnames=['e^entropy']);
Explanation: Finally, I run the MC algorithm over on the model using the algorithm and sample size I've observed in the learning materials.
End of explanation
<END_TASK>