Python: Sampling from a normal population

Here is code to create N samples of size n from a normal population with mean mu and standard deviation sigma.

"""
File             normalsim.py
Author           Dr. Ernesto P. Adorio
Desc             Sampling from the normal distribution.
Version          0.0.1 Sep. 1, 2011
"""



import numpy as np
import scipy.stats as stat
from math import sqrt


def  normalsim(mu, sigma, n, N, statistic="mean"):
    out = [0] * N
    for i in range(N):
        sample = stat.norm.rvs(mu, sigma, size=n)
        #print sample, 
        if statistic=="mean":
           out[i] = np.mean(sample)
        elif statistic=="svar":
           out[i] = np.var(sample, ddof=1)
        elif statistic=="pvar":
           out[i] = np.var(sample)
        elif statistic=="ssdev":
           out[i] = sqrt(np.var(sample, ddof=1))
        elif statistic=="psdev":
           out[i] = np.var(sample)
        elif statistic== "min":
           out[i] = min(sample)
        elif statistic == "max":
           out[i] = max(sample)
        else:
           raise ArgumentError,"bad statistic."
        #print out[i]
    return out

v = normalsim(10, 4, 10, 1000, statistic="svar")
print np.mean(v), np.var(v)

When I run the program under the time command, time python normalsim.py, the shell outputted

15.7770857263 55.2458134956

real    0m2.681s
user    0m1.080s
sys     0m0.070s

Problem: Can you rewrite the code so it will be faster? HINT: a comparison is always used inside the loop!

Python: converting names to surname, firstname format.

we want to convert names like

Jose P. Rizal to Rizal, Jose P.
Jose P. Rizal Jr. to Rizal Jr.,   Jose P.
Ma. Jose P. Rizal III to Rizal III, Ma. Jose P.

Here is Python code to perform the conversion. I am wondering if one can write a more elegant version, where by elegant, it does the conversion in a faster, shorter way!

"""
file     lastfirstname.py
author   Ernesto P. Adorio, PhD.
desc     Converts name to last name, first name format.
version  0.0.1 july 27, 2011
"""

def makestr(sarray, lastindex):
    """
   
    """
    output = ""
    for s in sarray[:lastindex]:
       output += (" " + s )
    return output

def surnamefirstname(name):
    """
    Translates a name in firstname surname to
          surname, firstname format.
    Complications arise because surname can include "Jr." or romanized numbers "III"
    """
    parts = name.split()
    n = len(parts)
    if n == 2:
       # easiest case
       return parts[1] + ", " + parts[0]
    elif n > 2:
       if parts[-1].lower() in ["jr.", "jr"]:
          return parts[-2] +" Jr.," + makestr(parts, -2)
       elif parts[-1].lower() in ["ii", "iii", "iv", "v"]:
          return parts[-2] + " " +parts[-1].upper()+ ","  + makestr(parts, -2)
       else:
          return parts[-1] + "," + makestr(parts, -1)
 
print surnamefirstname("Ma. Jose P. Rizal III")

 
       

Post you own better version in the comments. We shall learn together to improve our Python skills.

A beginning Python programming challenge: replace slow array.index() calls.

For simplicity, consider the following arrays:

i      0  1   2  3  4  5  6  7  8  9  10

order1 6  9   7 10  5  8 3  2  4   0   1

order2 7  6  10  8  5  9 1  4  2   0   3

Now consider the following Python snippet of code:

myindex = range(11)

for i in range(11):

    myindex[i] = order2[order1.index[i]]

As an example consider a current i value of 0. The value of index1[0] is 6. The position of 6 in order2 is of 6 in index2 array is 1(first element has index0). Therefore myindex[0] = 1.

On the other hand, the value of index1[1] is 9. The ordinal position of 9 in index2 is 5. Hence myindex[1] = 5.

Problem 1. Can you replace the code without using dictionaries and the index function?
You may not modify inplace the order1 and order2 arrays. You may create an additional extra array.

Problem 2.  Rewrite the code to use a dictionary. This is much easier.

Imagine your are doing the for loop of 1,000,000 array items! Have fun!

Python: Generating all subsets of size k for n objects.

In our main blog, we presented a routine for generating all subsets of set, expressed as an array or vector of indices to the n objects (0 to n-1).See Generating all subsets of a set with n objects In this article we present an "odometer" method of generating all subsets of specific size k where k < n. For example, for n = 5 and k = 3, we expect the function to return the set containing the following vectors: [0,1,2], [0, 1,3], [0,1,4], [0,2,3] and so on up to [2,3,4]. Here is our Python code.

# -*- coding: utf-8 -*-



"""

File      ksubsets.py

Author    Ernesto P. Adorio, Ph.D.

UPDEPP, UP Clarkfield  

Universty of the Philippines Extension Program in Pampanga

Desc      generates subsets (as indices) with k elements from a set of size n.

Version   0.0.1  2011.07.03

"""





def  ksubsets (n,  k):

"""

generates the list of subsets (as index array) of n elements

chosen k at a time. Combinations actually!

"""

#sanity check

if n <0 or k <0 or n - k < 0:
        return None
    N = range(k)
    retval = [N[:]]
    print "initial N=", N
    
    #update current N vector.
    level = k-1
    while level >=0:

d = N[level]

#test if element is already  at maximum value.

if d < n - k + level:  
            # No, increment the value.  
            N[level]+= 1
            #repair trailing elements.
            for i in range(level+1,  k):
                N[i] = N[i-1] + 1
            #append newly generated subset.
            retval.append(N[:])
            level = k-1
        else:
            # position N[level] already at maximum value!
            level -= 1
    return retval
    
for i,  subset in enumerate(ksubsets(5,  3)):
    print i,  subset

When the above program is run, it outputs the expected results:

0 [0, 1, 2]
1 [0, 1, 3]
2 [0, 1, 4]
3 [0, 2, 3]
4 [0, 2, 4]
5 [0, 3, 4]
6 [1, 2, 3]
7 [1, 2, 4]
8 [1, 3, 4]
9 [2, 3, 4]

How does the code work? First it initializes the first member of the set of subsets as
[0, 1, ..., k-1]. Next, starting from the rightmost (level = k-1)position it checks whether the current element is at maximum. If it is not it is incremented and any trailing positions adjusted to be 1 greater than at the preceding position. Otherwise the level is decremented and the process repeats.

Exercise: write a generator version of the above code. We will present the solution on the next revision.

The itertools module has a combinations function. After importing, one simply writes:
list(combinations (range(5), 3)).. Python just makes it TOO easy!

Python in Electrical Engineering: Converting between Y and Delta impedances

Given the branch impedances in delta configuration, what is the equivalent circuit impedance if connected in a Y configuration? the folowing short Python codes ease the conversion between Y and Delta impedances.

"""
file       deltay.py
author     Ernesto P. Adorio, Ph.D.
           UPDEPP (U.P. at Clarkfield)
desc       translate delta impedances to Y impedances.
version    0.0.1    04.21.2011
"""



def delta2yZ(Zab ,Zbc, Zca):
    """
    translate delta impedances to Y impedances.
    Arguments
       Zab, Zbc, Zca - complex impedances of Deta connection.
    Output
       Za, Zb, Zc - equivalent Y connection impedances to commom point.
    """
    denom = Zab + Zbc + Zca
    return (Zca * Zab/denom,Zab* Zbc/denom, Zbc * Zca/ denom)

def y2deltaZ(Za, Zb, Zc):
    """
    translate Y impedances to delta impedances.
    Arguments
       Za, Zb, Zc - complex Y connected impedances 
    Output
       Zab, Zbc, Zca - equivalent delta connection branch impedances.
    """
    num = Za*Zb + Zb*Zc + Zc* Za
    return (num/ Zc, num / Za, num/Zb)

We will continue this with an example later.

Exercise: State the formulas used by the above routines y2deltaZ(), delta2yZ() routines. It is so easy to write the functions that we did not bother to write the documentation for this first version, but we will later!

Python in Elecrical Engineering: Computing Symmetrical Components.

Given three-phase voltages V1, V2, V3 in phasor (r, phi) form, where r is the magnitude of the phasor, and phi is the phase angle. the following formulas compute the corresponding symmetrical components:

Null sequence component: $$E_0 = (1/3)(V1 + V2 + V3)$$.

Positive sequence component: $$E_1 = (1/3)(V1 + aV2 + a^2 V3)$$.

Negative sequence component: $$E_2 = (1/3)(V1 + a^2V2 + a V3)$$.

Here the a is the operator which rotates the phasor by 120 degrees. Power engineers usually specify the input voltages or current in phasor form. The formula shows we have to add phasors termwise and thus the routine to compute symmetrical components first convert them to complex numbers. Complex numbers are already handled naturally in Python! and here is the full source code:

"""
file    symcomp.py
author  Ernesto P. Adorio, Ph.D.
UPDEPP (U.P. Clarkfield)
desc    basic symmetrical components for three phase phasors.
Phasors are tuples of the form (r, phi) where r and phi are the
magnitude(absolute value) and phi is the phase angle.
version 0.0.1 april 20, 2011
"""


from math  import *
from cmath import *

zTOL = 1.0e-8
DTOR = pi / 180.0
RTOD = 180.0 / pi



ONETWENTYRAD = 120 * DTOR
_a_  = rect(1, 120 * pi/ 180.0)
_a2_ = _a_.conjugate()

def dtor(degree):
"""
Converts degree to radians.
"""
return degree * DTOR

def rtod(radian):
"""
Converts radians to degree.
"""
return radian * RTOD

def  z2pair(z):
"""
returns z in pair form.
""" 
return (z.real, z.imag)

def  pair2z(re, im):
"""
returns a complex number from the components. 
"""
return re + im*1j

def  polar2pair(v):
"""
(r, phi) to (re, im) pair form.
"""
# Special checking for zero imaginaries!
# some computations result in an additional 0j
r, phi = v
if type(r) == type(1j):
r =r.real
if type(phi) == type(1j):
phi = phi.real

z = rect(r, phi)
return z.real, z.imag

def  pair2polar(re, im):
return polar(pair2z(re,im))

def  zround(v):
if abs(v[0]) < zTOL:
        return (0.0, 0.0)

def  a(v):
     """
     Applies the a operator to a phasor v in polar form.
     It adds a 120 degree to the phase of the phasor v.
     if v is real, the angle is zero.
     """
     if type(v) != type((0,0)):
        v = (v, 0)
     r     = abs(v[0])
     theta = v[1]

     newangle = theta + ONETWENTYRAD
     return (r * cos(newangle), r * sin(newangle))

def a2(v):
     """
     Applies the a operator to a phasor v in polar form.
     It adds a 120 degree to the phase of the phasor v.
     """
     if type(v) != type((0,0)):
        v = (v, 0)
     r     = abs(v[0])
     theta = v[1]

     newangle = theta + 2*ONETWENTYRAD
     return (r * cos(newangle), r * sin(newangle))


def  symcomp(v1, v2, v3):
     """
     Returns the symmetrical components of the three phasors v1, v2, v3
     which are in tuple (r, theta) form.
     """
     #Convert first to complex rectangular form.
     v1z  = rect(v1[0], v1[1])
     v2z  = rect(v2[0], v2[1])
     v3z  = rect(v3[0], v3[1]) 

     av2  = _a_ * v2z
     a2v2 = _a2_ * v2z
     av3  = _a_* v3z
     a2v3 = _a2_* v3z

     #Null sequence  component.
     E0 = polar((v1z.real+ v2z.real+  v3z.real)/3.0 +  (v1z.imag+ v2z.imag+ v3z.imag)/3.0*1j)

     #Positive sequence component.
     E1 = polar((v1z.real + av2.real +  a2v3.real)/3.0+ (v1z.imag+ av2.imag+ a2v3.imag)/3.0*1j)

     #Negative sequence component.
     E2 = polar((v1z.real + a2v2.real +  av3.real)/3.0+ (v1z.imag+ a2v2.imag+ av3.imag)/3.0*1j)

     return (E0, E1, E2)


def  symcompz(v1, v2, v3):
     """
     Returns the symmetrical components of the three phasors v1, v2, v3
     which are in tuple (r, theta) form.
     """
     av2  = _a_ * v2
     a2v2 = _a2_ * v2
     av3  = _a_* v3
     a2v3 = _a2_* v3

     #Null sequence  component.
     E0 = polar((v1.real+ v2.real+  v3.real)/3.0 +  (v1.imag+ v2.imag+ v3.imag)/3.0*1j)

     #Positive sequence component.
     E1 = polar((v1.real + av2.real +  a2v3.real)/3.0+ (v1.imag+ av2.imag+ a2v3.imag)/3.0*1j)

     #Negative sequence component.
     E2 = polar((v1.real + a2v2.real +  av3.real)/3.0+ (v1.imag+ a2v2.imag+ av3.imag)/3.0*1j)

     return (E0, E1, E2)

def  symcomp2phasors(E0, E1, E2):
     """
     Recreates the phasors form the symmetrical components.
     """
     V1 = polar(rect(E0[0], E0[1]) + rect(E1[0], E1[1]) + rect(E2[0], E2[1]))
     V2 = polar(rect(E0[0], E0[1]) + _a2_* rect(E1[0], E1[1]) + _a_ *rect(E2[0], E2[1]))
     V3 = polar(rect(E0[0], E0[1]) + _a_* rect(E1[0], E1[1]) + _a2_ * rect(E2[0], E2[1]))
     return V1, V2, V3


if __name__ == "__main__":
   #extreme cases.
   I1 = polar(10)
   I2 = polar(0)
   I3 = polar(0)
   print symcomp(I1, I2, I3)   
  
  

   #extreme cases, balanced system.
   i1 = 1
   i2 = -0.5+sqrt(3)/2.0 * 1j
   i3 = -0.5-sqrt(3)/2.0 * 1j

   I1 = polar(i1) 
   I2 = polar(i2)
   I3 = polar(i3)
   E0, E1, E2 = symcomp(I1, I2, I3)   
   print "original phasors=", I1, I2, I3
   print "symmetrical components:", E0, E1, E2
   phasors = symcomp2phasors(E0, E1, E2)
   print "recovered phasors:", phasors
   #include more here! from published books or other sources.

We have not yet fully tested the above code, and we would appreciate it if our readers point out any any errors in the program. When the current version 0.0.1 is run, it outputs the following:

$ python symcomp.py
((3.3333333333333335, 0.0), (3.3333333333333335, 0.0), (3.3333333333333335, 0.0))
original phasors= (1.0, 0.0) (0.99999999999999989, 2.0943951023931957) (0.99999999999999989, -2.0943951023931957)
symmetrical components: (7.4014868308343765e-17, 3.1415926535897931) (1.295260195396016e-16, 0.0) (0.99999999999999989, 0.0)
recovered phasors: ((1.0, 9.0639078007724287e-33), (0.99999999999999978, 2.0943951023931953), (0.99999999999999978, -2.0943951023931953))
$

Note that in the second example, the null-sequnce component is the zero vector (zero magnitude). Out routines are able to recover the original phasors from the computed symmetrical componsnts.

Using javascript for data visualization and animations.

Visit neurosoftware.ro.


We are still learning how to run jsxscript in our Wordpress blog.