# We need to work around not being allowed to have a negative exponent
# using scipy's powerlaw package. Thank you stackoverflow and Jeff Alstott
#https://stackoverflow.com/questions/17882907/python-scipy-stats-powerlaw-negative-exponent
# Here is Alstott's Python package powerlaw that does this
#https://github.com/jeffalstott/powerlaw
!pip install powerlaw
import powerlaw
import matplotlib.pyplot as plt
import numpy as np
a, xmin = 3.0, 1.0 #Power law starts at xmin = 1 with slope of -3
N = 10000
vrs = powerlaw.Power_Law(xmin = xmin, parameters = [a],).generate_random(N)
# plotting the PDF estimated from variates
bin_min, bin_max = np.min(vrs), np.max(vrs)
bins = 10**(np.linspace(np.log10(bin_min), np.log10(bin_max), 100))
counts, edges = np.histogram(vrs, bins, density=True)
centers = (edges[1:] + edges[:-1])/2.
# plotting the expected PDF
xs = np.linspace(bin_min, bin_max, 100000)
plt.plot(xs, [(a-1)*xmin**(a-1)*x**(-a) for x in xs], color='red')
plt.plot(centers, counts, '.')
plt.xscale('log')
plt.yscale('log')
plt.savefig('powerlaw_variates.png')

```
Requirement already satisfied: powerlaw in /root/venv/lib/python3.7/site-packages (1.4.6)
Requirement already satisfied: matplotlib in /shared-libs/python3.7/py/lib/python3.7/site-packages (from powerlaw) (3.4.2)
Requirement already satisfied: scipy in /shared-libs/python3.7/py/lib/python3.7/site-packages (from powerlaw) (1.7.0)
Requirement already satisfied: numpy in /shared-libs/python3.7/py/lib/python3.7/site-packages (from powerlaw) (1.19.5)
Requirement already satisfied: mpmath in /shared-libs/python3.7/py/lib/python3.7/site-packages (from powerlaw) (1.2.1)
Requirement already satisfied: pyparsing>=2.2.1 in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from matplotlib->powerlaw) (2.4.7)
Requirement already satisfied: cycler>=0.10 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from matplotlib->powerlaw) (0.10.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from matplotlib->powerlaw) (1.3.1)
Requirement already satisfied: pillow>=6.2.0 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from matplotlib->powerlaw) (8.3.1)
Requirement already satisfied: python-dateutil>=2.7 in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from matplotlib->powerlaw) (2.8.1)
Requirement already satisfied: six in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from cycler>=0.10->matplotlib->powerlaw) (1.16.0)
WARNING: You are using pip version 21.1.2; however, version 21.1.3 is available.
You should consider upgrading via the '/root/venv/bin/python -m pip install --upgrade pip' command.
```

#Now create a uniform distribution
from scipy.stats import uniform
vrs0 = uniform.rvs(size=10000)
plt.hist(vrs0, density=True, histtype='stepfilled', alpha=0.2)
plt.legend(loc='best', frameon=False)
plt.show()

```
No handles with labels found to put in legend.
```

#Now plot power law vrs
plt.hist(vrs, density=True, histtype='stepfilled', alpha=0.2)
plt.legend(loc='best', frameon=False)
plt.show()

```
No handles with labels found to put in legend.
```

#Oops need those nice log axes for a power law
plt.hist(vrs, density=True, bins = np.logspace(start = np.log10(1), stop = np.log10(20)), histtype='stepfilled', alpha=0.2)
plt.legend(loc='best', frameon=False)
plt.gca().set_xscale("log")
plt.show()

```
No handles with labels found to put in legend.
```

#Power law scaling on both axes please?
plt.hist(vrs, density=True, bins = np.logspace(start = np.log10(1), stop = np.log10(20)), histtype='stepfilled', alpha=0.2)
plt.legend(loc='best', frameon=False)
plt.gca().set_xscale("log")
plt.gca().set_yscale("log")
plt.show()

```
No handles with labels found to put in legend.
```

#Just show these two sets of random variables
# The uniform ones are showing but they are down between 0 and 1.
plt.plot(vrs0)
plt.plot(vrs)
plt.show()

#Do better ... use subplots
fig,axs = plt.subplots(2)
axs[0].plot(vrs0)
axs[1].plot(vrs, 'tab:orange')

#Combine these numbers a different way so two distributions form one long vector
vrscatted = np.concatenate([vrs0, vrs])

#Now plot distribution
plt.hist(vrscatted, density=True, bins = np.logspace(start = np.log10(.01), stop = np.log10(80)), histtype='stepfilled', alpha=0.2)
plt.gca().set_xscale("log")
plt.show()

# This will write file. I have turned off row numbers from being written with "index=False"
# So one can strip them off later?
import pandas as pd
df = pd.DataFrame(vrscatted)
df.to_csv('vrscatted.csv', index = False)

print(vrscatted)

```
[0.38458859 0.70233379 0.61619905 ... 1.6368232 1.42320643 2.34780992]
```

min(vrs)
max(vrs0)

len(vrs)

from numpy import genfromtxt
temporary = genfromtxt('vrscatted.csv')

from matplotlib.pyplot import plot
plot(temporary)