Skewness and Kurtosis

http://www.dreamincode.net/code/snippet3048.htm

 

Found a code for skewness and kurtosis. Just putting it here for reference.

from math import *
from Numeric import *
n=input("Enter the number of class:")
h=input("Enter class interval:")
sum_f=0
sum_fm=0
total1=0
total2=0
total3=0
a=zeros(n,Int)
b=zeros(n,Int)
for i in range(n):
    print "For class-%s"%(i+1)
    f=input("Frequency:")
    l=input("Lower limit:")
    u=l+h    #upper limit
    ul=(l+u)/2.0	#mean of upper and lower limit
    ful=f*ul	#multiplication of frequency and mean
    a[i]=ul
    b[i]=f
    sum_f=f+sum_f
    sum_fm=ful+sum_fm
x=sum_fm/sum_f		#arithmetic mean
for i in range(n):
    d=a[i]-x
    v1=pow(d,2)*b[i]
    v2=pow(d,3)*b[i]
    v3=pow(d,4)*b[i]
    total1=v1+total1
    total2=v1+total2
    total3=v1+total3    m2=total1/sum_f    #2nd moment
m3=total2/sum_f    #3rd moment
m4=total3/sum_f    #4th moment
s=pow(m3,2)/pow(m2,3)    #skewness
k=m4/pow(m2,2)    #kurtosis
print "The skewness is %s & kurtosis is %s"%(s,k)

Interpreting FFT results and plotting

I don’t frequently use FFT, so it still confounds every time I try to interpret the results from FFT, I’ll try to organize my understanding here…

INFORMATION ABOUT THE DATA FILE TO PROCESS

data = array of samples
data_pts = total samples in data file (number of samples/lines in data file)
data_srate = sampling rate (points per second)
data_stime = total time covered in seconds (data_pts / data_srate)

THE numpy.fft() FUNCTION AND THE RESULTS

fft_result = fft(data)
fft_freqs = fftfreq( len(fft_result) )

(fft_result has a real and imag part, but I’m focusing on the real part here)

len(fft_result) == len(data)

but this is mostly meaningless

fft_result[n] corresponds to fft_freqs[n]

PRECISION
The interval difference if fft_freqs equals the inverse of data_stime.
This interval has nothing to do with the number of samples which is what confused me most. I had expected that having 1024 samples per second would return a more precise frequency plotting versus 2 samples per second.
In other words, if the sampled time period is 10 seconds long. No matter if you have 2 or 100 samples per second, the intervals of the frequency ‘bins’ will be 0.1. The number of 0.1-bins depends on the number of samples. In this case, 2 samples per second will produce 20 bins ranging from -1.0 to 0.9; 100 samples will give 1000 bins, from -50.0 to 49.9.

dfreq = 1/data_stime
range = (-dfreq*data_stime*data_srate/2 , dfreq*data_stime*data_srate/2 - dfreq)
      = (-data_srate/2                  , data_srate/2 - dfreq                 )

IN SHORT
That seems more complicated than it needs to be, but it was the thought process that helped me arrive at a better solution. After doing the math and comparisons I find that there is no need to do all these calculations. The sampling rate only needs to be applied to the freq bins, similar to dividing the x-axis array by the rate when plotting the data. Squaring the FFT values helps accentuate the peaks.

    data = []       #array to hold data read from file
    srate=2.        #sample rate
    samples = 0     #counter for the number of lines in file (# of data samples)
    with open(filename, 'r') as rawfile:
        for i, line in enumerate(rawfile):
            data.append(float(line))
        samples = i + 1

    x = arange(samples)/srate

    # PLOT THE ORIGINAL DATA READ FROM FILE
    plt.plot(x, data)
    plt.show()

    # APPLY FFT
    FFT = fft(data)
    freqs = fftfreq(samples) * srate
    print 'bin range=', [freqs.min(), freqs.max()]
    print 'bin interval=', srate/samples

    # PLOT THE POSITIVE HALF OF FREQ VALUES
    plt.plot(freqs[:samples/2], (abs(FFT.real)**2)[:samples/2], 'g*')
    plt.show()

Ax = B is not the same as x = A(inv) B

When given the task to change my explicit ocean model to an implicit version, I had assumed that I would need to invert the ‘A’ matrix to solve Ax=B, and then do some matrix math. I found a few web sites that warn against that because it is possibly a lot more computationally intensive than solving directly. John Cook’s blog explains more about why.

Python’s NumPy has linalg.solve(A, B), which returns the ‘x’ array

x = numpy.linalg.solve(A,B)

It uses a LU decomposition method for solving (not inversion).

Sliding plates implicit

The next assignment was to rewrite the sliding plates using an implicit scheme. Here is a video of the animation produced using Visual Python arrows.

Colored according to speed. The initial velocity of the fluid is zero and the upper plate moves at a constant speed.
The visualization setup using Visual Python looks like this:

#SET UP VISUALIZATION SCENE
scene = display(title='Velocity profile', width=800, height=800, center=(u/2,y_divs/2,0), 
                autoscale=True, fullscreen=False)
base = box(pos=(u/2,-0.5,-2), length=u, height=0.4, width=4)
plate = box(pos=(u/2,y_divs+1.5,-2), length=u, height=0.4, width=4)
quiver = [] #TO HOLD ALL THE VELOCITY VECTOR ARROWS
quiver.append( arrow(pos=(0,0,0), axis=(BOUNDS[0,0],0,0),
                         shaftwidth=0.5, color=color.blue) )
for each in xrange(y_divs):
    quiver.append( arrow(pos=(0,each+1,0), axis=(UM[each,0],0,0),
                         shaftwidth=0.5, color=color.blue) )
quiver.append( arrow(pos=(0,y_divs+1,0), axis=(BOUNDS[-1,0],0,0),
                         shaftwidth=0.5, color=color.red) )
scene.autoscale = False

The main loop with the visualization update looks like this:

    for t in arange(t_divs-1): #START FROM 1, SKIP ZERO
        B = zeros(y_divs)
        B[:] = -1 * UM[:,t]
        B[ 0] -= BOUNDS[ 0,t+1] * alpha
        B[-1] -= BOUNDS[-1,t+1] * alpha
        
        UMnext = linalg.solve( A, B )
        UM[:,t+1] = UMnext[:]
        
        for y in arange(y_divs):
            #UPDATE VECTOR LENGTH BASED ON SPEED
            quiver[y+1].axis = (UM[y,t+1],0,0)
            #CHANGE ARROW COLOR BASED ON SPEED: Blue is slow, Red is fast
            quiver[y+1].color = (UM[y,t+1]/u,0.3,1.-UM[y,t+1]/u)

Sliding Plates

Completed the sliding plates assignment.
Assignment: Simulate the vertical velocity profile of a liquid between two infinitely long plates separated by 40mm. The bottom plate remains fixed, but the top plates begins moving at a constant 40m/s in the +x direction. Using the explicit finite difference method.

Image

Visual python makes it so easy to make a nice graphical simulation. I use arrows along the profile that change from blue to red as they increase in length to show the velocity.

Image

tidbits from reading

.index() vs. where()

.index() is the builtin function for finding the index of something in a list or tuple. This does not work with NumPy arrays! NP arrays rely on where() to return an array of indices that match the condition.

—————————
for-else statement

‘else’ can be used after a for-loop. But here, the ‘else’ statement only executes if the for-loop completes all iterations. Think of it as using the for-loop to search for something, or expecting a break from the loop. If nothing is resolved or no break occurs then the else statement can be used for notification. The else can be used after a try-except in a similar way; executing only when nothing breaks down, throwing an exception.

—————————
how to prevent a new-line character in a print() statement

Use a comma at the end. You can have consecutive print() statements print to the same line simply by adding a comma at the end of each line.

—————————
Use a dictionary for formatting a long output string

%20s will format a string in the list after a text. But when using this formatting method, the order must match in the supplied list. %(keyword)20s or %(keyword)10f with a supplied dictionary can avoid the ordering issue and having to have the same variable or value in a list multiple times.

—————————
use an asterisk * in front of a functions input argument to accept an infinite number of args.

def MyFunction( *args ): This will accept any number of arguments and create a list to use within the function.

—————————
ternary operation or inline if (iif)

In Java you can use a iif like this ( x == 2 ? y = 4 : y = 0 ).
( condition ? if true do this : if false do this)
I found that Python also has a version of this but ordered differently:
(do this) if (condition) else (do this if false)
Not the same order as java but still pretty straight forward.

print(“Hello world!”)

For the first entry of my programming activities record , I will start by posting my current script template. I plan on periodically updating this as I learn more about Python, but I guessing it won’t change much from this format. Maybe I should put it on a page instead of post.

_____________template.py_______________

#!/usr/bin/env python
# -*- coding: utf-8 -*-
“””
Created on %(date)s

@author: Ripley6811

@abstract:

@input: none

@output: none

@acknowledgments:
“””
#from numpy import * #imports ndarray(), arange(), zeros(), ones(), and other array methods
#from visual import * #imports NumPy.*, SciPy.*, and Visual objects (sphere, box, etc.)
#import matplotlib.pyplot as plot #use plot(x,y) and show() for plotting
#from pylab import * #imports NumPy.*, SciPy.*, and matplotlib.*

def main():

“””Description of main()”””

if __name__ == ‘__main__’:

main()

”’######################################
linspace(0,2,9) > array of 9 numbers from 0 to 2
zeros( (10,10,10) ) > a 3-d array of 0.0 (pass a list)
ones( (3,4,5) ) > a 3-d array of 1.0’s
”’######################################