Power Law Distributions
=======================

(Note: Be sure to either download the data that goes with this notebook or uncomment/edit the lines below needed to get your own data before executing the code in this notebook.)

It turns out that "power law" distributions are ubitquitous in nature[1]. 

1 <http://arxiv.org/abs/cond-mat/0412004>

In the cells that follow I measure the number of files of different sizes on my hard disk. Then I proceed to quantify (in the simples possible way) the distribution of file sizes in a certain range as a power law distribution.



```python
%pylab inline

import pandas as pd                    # pandas is the "python data analysis package"
import os                              # operating system services
from scipy.optimize import curve_fit   # non-linear curve fitting
```

    Populating the interactive namespace from numpy and matplotlib



```python
def getFileSizes(folder, minsize=2000, maxsize=10**9):
    """
    return a list of file sizes in a folder
    """
    sizes = []
    for (path, dirs, files) in os.walk(folder):
        for file in files:
            try:
                filename = os.path.join(path, file)
                size=os.path.getsize(filename)
                if size>minsize and size<maxsize:
                    sizes.append(size)
            except OSError:
                pass
    return sizes
```


```python
if 0: # set this to "1" to collect data
    minsize=10**4
    maxsize=10**10
    sizes = getFileSizes('/Users/steve/Desktop',minsize,maxsize) # <--- edit this to match some folder with a lot of files
    print(len(sizes), 'files between', minsize, 'and', maxsize)
```


```python

if 0: # set this to "1" to write the file
    df=pd.DataFrame({'file_sizes':sizes})
    df.to_csv('myData.csv', sep=',')
else:
    df=pd.read_csv('myData.csv')
    df.head()

```


```python
ns,bins,patches = hist(log(array(sizes)+1),40)
```


```python
xvals=(bins[:-1]+bins[1:])/2.0
counts=ns+0.5
yvals=log(counts)
title("distribution of file sizes")
ylabel("log(count)")
xlabel("log(size)")
plot(xvals,yvals,'b.')
```

Fitting Data with Matrices 
===========================

We can see that we have something that looks like it might be a power law. How can we tell? We'll use "linear least squares" analysis.

When fitting some experimental data to a model, there are typically
a set of paramters $(\alpha_1, \alpha_2, \ldots, \alpha_n)$ that
need to be found so that the function $\chi^2$ is a minimum. 
We can think of the paramters as an n-vector called $\alpha$. The
data consist of a set of paired numbers $(x_i, y_i)$ and we can
define $\chi^2$ as 

\begin{equation}
\chi^2 = \sum_{i=1}^{m} {(y_i - y^\star(x_i,\alpha))^2\over \sigma_i^2},
\end{equation}

where $\sigma_i$ is the uncertainty in $y_i$, and $y^\star$ is the
model value corresponding to $x_i$ (and the value of the parameters
$\alpha$). If the model $y^\star(x,\vec\alpha)$ is *linear* then we
can write $y^\star$ as:

\begin{equation}
y^\star (x) = \alpha_1 X_1(x)+\alpha_2 X_2(x)+\ldots+\alpha_n X_n(x).
\end{equation}

or

\begin{equation}
y^\star(x) = \sum_{j=1}^{n} \alpha_j X_j(x),
\end{equation}

where the $X_j$ are a set of functions that are to be superposed
with coefficients $\alpha_j$ to form the model. A common example
of this would be a straight line fit $y^\star(x) = mx+b$ where the
functions are $x$ and 1 for the paramters $m$ and $b$ respectively.

Notice that, since the $x$'s are at fixed measured values, this whole thing can be 
cast into a single matrix equation. If we define the following $M$ matrix:

\begin{equation}
M_{ij} = X_j(x_i)
\end{equation}

and if we think of $y^\star$ as an $m$-vector of model 
values evaluated at the $m$ data points $x_i$ and $y$ as an 
experimental $m$-vector $y$, then the model vector
can be evaluated by matrix multiplication as:

\begin{equation}
y^\star = M \alpha
\end{equation}

and the difference vector between the model and the data becomes simply

\begin{equation}
y - y^\star = y - M\alpha
\end{equation}

Finally we can write $\chi^2$ as 

\begin{equation}
\chi^2 = (y - M\alpha)^T S (y - M\alpha) \,\,\,\,\,\,\,\,(1)
\end{equation}

where $S$ is a diagonal matrix of $1/\sigma_i^2$ and
the $T$ on the left side means "transpose".

Now it is a simple matter to find the values of the paramters that
minimize $\chi^2$. We require that $\nabla_{\alpha}\chi^2 = 0$. 
If you play around with this you can convince yourself that this means:

\begin{equation}
M^T S(y - M\alpha) = 0
\end{equation}

and therefore:

\begin{equation}
\alpha = (M^T SM)^{-1}M^T S y.  \,\,\,\,\,\,\,\,(2)
\end{equation}


```python
#
# define the function values (these are the columns of the "M" matrix
#
# generate arrays that contain "functions" of x that are combined to form the model
#

def doFit(funcs, xvals, yvals, sigma):
    """
    doFit uses matrices to solve for the linear least squares 
    fit parameters based on the list of 'funcs', and the x-y values
    and their corresponding uncertainties (sigma).
    
    returns the array of model parameters
    """
    #
    # define the S matrix in terms of uncertainties
    #
    
    S=diag((1.0/sigma)**2)
    
    #
    # define the model matrix
    #
    
    M = array(funcs).T
    MtM = M.T.dot(S.dot(M))
    MtMInv = inv(MtM)
    MtY = M.T.dot(S.dot(yvals))
    alpha = MtMInv.dot(MtY)      # see eq 2 above.
    
    #
    # solve for the parameters, and return
    #
    ystar = M.dot(alpha)
    
    return (alpha, ystar, MtMInv)

fx1 = ones(len(xvals))   # X_1(x) = 1
fx2 = xvals              # X_2(x) = x

sigma=1.0/sqrt(abs(counts))

alpha, ystar, fcov = doFit([fx1, fx2], xvals, yvals, sigma)

mFit = alpha[1]
bFit = alpha[0]

errorbar(xvals, yvals, sigma, fmt='r.') 
title("File Sizes on My HD")
ylabel("log(count)")
xlabel("log(size)")
plot(xvals,ystar, 'g-')

```

Uncertainties in the Parameters
-------------------------------

There are a variety of techniques for estimating the uncertainty in the parameters you compute. The approach we'll take is one of the more "brute force" approaches, but it's pretty robust. We're going to fabricate data that has the same statistics as our original data (do the degree that we know and/or care). Then we'll fit the fabricated data many times (1000 times in this case) and we'll look at the statistics of the fit parameters produced. From this we can tell how well our fit has done at narrowing down the parameters we computed.



```python
#
# Now do monte-carlo data fabrication and fit analysis
#

N=1000    # number of samples to take
mList=[]  # keep track of monte-carlo fit parameters
bList=[]

for i in range(N):
    """
    Generate mc data with the same statistical properties as the real data.
    Repeat the fit for each set, and record the parameters.
    """
    mcY = ystar + sigma*normal(size=len(xvals))
    mcAlpha,mcYstar,mccov = doFit([fx1,fx2],xvals,mcY,sigma)   # repeatedly fit mc data
    mList.append(mcAlpha[1])
    bList.append(mcAlpha[0])
    
#
# Compute the statistics of the mc-results
#

marr=array(mList)
barr=array(bList)

mAvg = marr.sum()/N
bAvg = barr.sum()/N
delM = marr-mAvg
delB = barr-bAvg
sigM = sqrt((delM*delM).sum()/(N-1))  # sigM is the std-deviation of the m values
sigB = sqrt((delB*delB).sum()/(N-1))  # sigB is the std-deviation of the b values

#
# plot the fit
#
print("Slope=", mFit, '+/-', sigM, "(", mFit - sigM,",",mFit + sigM, ")")
print("Intercept=", bFit, '+/-', sigB, "(", bFit - sigB,",",bFit + sigB, ")")

#
# Just for fun print out the diagonals of the covariance matrix.
#
print # print a blank link
print("Compare to cov-matrix for fun:")
print("sqsrt(cov[1,1]) (should be sigma m)", sqrt(fcov[1,1]))
print("sqrt(cov[0,0]) (should be sigma b)", sqrt(fcov[0,0]))

```

    Slope= -0.401012972104 +/- 0.00676272823764 ( -0.407775700342 , -0.394250243866 )
    Intercept= 10.2511101646 +/- 0.078149615394 ( 10.1729605492 , 10.32925978 )
    Compare to cov-matrix for fun:
    sqsrt(cov[1,1]) (should be sigma m) 0.00691993760803
    sqrt(cov[0,0]) (should be sigma b) 0.080079083965



```python
hist(marr) # look at the variation in "m" values
print("m average:", mAvg)
print("m sigma:",sigM)
```


```python
hist(barr)  # look at the variation in "b" values
print("b average:", bAvg)
print("b sigma:",sigB)
```


```python
#
# Just to see same data fit using the non-linear scipy package "curve_fit":
#

def fLinear(x, m, b):
    return m*x + b

popt, pcov = curve_fit(fLinear, xvals, yvals, p0=(alpha[1],alpha[0]), sigma=sigma)

m=popt[0]          # slope
dm=sqrt(pcov[0,0]) # sqrt(variance(slope))
b=popt[1]          # int
db=sqrt(pcov[1,1]) # sqrt(variance(int))
ystar=fLinear(xvals, m, b)

N=1000    # number of samples to take
mList=[]  # keep track of monte-carlo fit parameters
bList=[]

for i in range(N):
    """
    Generate mc data with the same statistical properties as the real data.
    Repeat the fit for each set, and record the parameters.
    """
    mcY = ystar + sigma*normal(size=len(xvals))  # generate fabricated data to fit
    mcpopt, mcpcov = curve_fit(fLinear, xvals, mcY, p0=(m,b), sigma=sigma)
    mList.append(mcpopt[0])  # store the fit paramters for the fab-data
    bList.append(mcpopt[1])
    
#
# Compute the statistics of the mc-results
#
marr=array(mList)
barr=array(bList)

mAvg = marr.sum()/N
bAvg = barr.sum()/N
delM = marr-mAvg
delB = barr-bAvg
sigM = sqrt((delM*delM).sum()/(N-1))  # sigM is the std-deviation of the m values
sigB = sqrt((delB*delB).sum()/(N-1))  # sigB is the std-deviation of the b values

errorbar(xvals, yvals, sigma, fmt='r.') 
title("File Sizes on My HD")
ylabel("log(count)")
xlabel("log(size)")
plot(xvals,ystar, 'g-')
print("Slope=", m, '+/-', sigM, "(", m-sigM,",",m+sigM, ")")
print("Intercept=", b, '+/-', sigB, "(", b-db,",",b+sigB, ")")

print() # print a blank link
print("Compare to cov-matrix for fun:")
print("sqsrt(pcov[0,0]) (should be sigma m)", sqrt(pcov[0,0]))
print("sqrt(pcov[1,1]) (should be sigma b)", sqrt(pcov[1,1]))

```

Project 7
=========

See this reference for other systems that exhibit power law behavior:

<http://arxiv.org/abs/cond-mat/0412004>

Or test your intuition and see if it's right. Many system exhibit this behavior!

Get some data, and test your idea. Use the matrix technique to do the fitting.


```python

```

ac0040da44bdffba94b688869e0040373b1dd3f1
ipynb
Jupyter Notebook
P07-Fitting-With-Matrices.ipynb
TejasAvinashShetty/sci-comp-notebooks
8ac7887a6c5da6fff8173febbc9a9968a81836af
MIT
P07-Fitting-With-Matrices.ipynb
TejasAvinashShetty/sci-comp-notebooks
8ac7887a6c5da6fff8173febbc9a9968a81836af
MIT
2021-02-02T23:14:34.000Z
2021-03-20T05:46:23.000Z
P07-Fitting-With-Matrices.ipynb
TejasAvinashShetty/sci-comp-notebooks
8ac7887a6c5da6fff8173febbc9a9968a81836af
MIT
Qwen/Qwen-72B
1. YES
2. YES