Handwritten number recognition with Python (Machine Learning)

Here I am again with Machine Learning! This time I’ve achieved a great result though (for me at least!). By using another great dataset from UCI I was able to write a decent ML script which scored 95% in the testing part! I am really satisfied with the result.

Here is a sample of what the script should be able to read (in the example the number 9):

Some numbers, as the one above, were clear, others not so clear, since they were handwritten and then somehow (I do not know how) converted into digital images.

I had a hard time figuring out how the attributes in the dataset were coded but in the end I managed to figure it out! I guess making up such a dataset was a really long and boring work.

Anyway here is my script and below you can find the result of the test on the last 50 numbers or so.

	import matplotlib.pyplot as plt
	import numpy as np
	import os
	# Import of support vector machine (svm)
	from sklearn import svm
	"""--------------------------------SETTINGS---------------------------------"""
	# Load data from .txt file
	os.chdir("C:\\")
	file = open("semeion.txt","r")
	data = file.read()
	file.close()
	# Get index max (number of rows in the data)
	def getIndexMax(data):
	dataSplitted = data.split("\n")
	return len(dataSplitted)
	# Prepare data for fitting function
	# Answer: if True, returns answers (labels)
	# Training: if True, returns training samples (with no answers, only samples)
	# Last: if True, returns only last training sample/answer. Useful for testing
	def returnDataToUse(data,index,answers=False,training=False,last=False):
	dataSplitted = data.split("\n")
	# Check that the index is not bigger than our dataset
	if index > len(dataSplitted):
	print("Index out of bounds, index max:",len(dataSplitted))
	return 0
	# This bit of code returns answers
	if answers and not training:
	firstLine = dataSplitted[0][:-1].split(" ")
	firstAnsw = firstLine[256:]
	firstAnsw2 = [int(i) for i in firstAnsw]
	firstAnswInt = firstAnsw2.index(1)
	correctAnswers = [firstAnswInt]
	i = 1
	while i < index:
	temp = dataSplitted[i][:-1].split(" ")
	temp2 = temp[256:]
	temp3 = [int(k) for k in temp2]
	temp4 = temp3.index(1)
	correctAnswers.append(temp4)
	i += 1
	completeAnswers = np.array(correctAnswers)
	if last:
	return completeAnswers[-1]
	else:
	return completeAnswers
	# This bit of code returns pure samples
	if training and not answers:
	firstLine = dataSplitted[0][:-1].split(" ")
	firstTraining = firstLine[:256]
	trainingArray = np.array([float(i) for i in firstTraining])
	i = 1
	while i < index:
	temp = dataSplitted[i][:-1].split(" ")
	temp2 = temp[:256]
	temp3 = np.array([float(k) for k in temp2])
	trainingArray = np.vstack((trainingArray,temp3))
	i += 1
	if last:
	return trainingArray[-1]
	else:
	return trainingArray
	# This function displays the image of the number (sample at row x)
	# and prints the answer the predictor should give us back
	def displayImage(data,row):
	# Split each row
	dataSplitted = data.split("\n")
	# Get the 'rowth' row
	strings = dataSplitted[row]
	# Split row into numbers(string), and avoid blank at the end
	stringsSplitted = (strings[:-1]).split(" ")
	# Get target and convert it into numbers, then in a numpy array
	risp = stringsSplitted[256:]
	rispInt = [int(i) for i in risp]
	rispNp = np.array(rispInt)
	# Print original data and number to guess in readable format
	print(rispInt)
	print("Number to guess:",rispInt.index(1),"\n")
	# Training array converted into float numbers
	training = stringsSplitted[:256]
	trainingFloat = [float(i) for i in training]
	# Building 16x16 (image) array
	#.
	i = 16
	k = 0
	img = np.array(trainingFloat[:16])
	while i <= len(trainingFloat):
	#print(i)
	#print(k)
	temp = np.array(trainingFloat[k:i])
	img = np.vstack((img,temp))
	k = i
	i += 16
	# Plot image
	plt.imshow(img,cmap=plt.cm.gray_r,interpolation="nearest")
	plt.show()
	"""------------------------------TRAINING------------------------------------"""
	# FIX THE NUMBER OF TRAINING SAMPLES
	trainingSamples = 1500
	# Gamma: gradient descent parameter
	clf = svm.SVC(gamma=0.01, C=100)
	# Index max
	print("Maximum index:",getIndexMax(data),"\n")
	answerArray = returnDataToUse(data,trainingSamples,answers=True)
	trainingAr = returnDataToUse(data,trainingSamples,training=True)
	x,y = trainingAr,answerArray
	#Fit the data
	print("Training...")
	clf.fit(x,y)
	"""------------------------------Sample prediction--------------------------"""
	# CHOOSE AN EXAMPLE TO PREDICT
	example = 1555
	predictQ = returnDataToUse(data,example,training=True,last=True)
	predictA = returnDataToUse(data,example,answers=True,last=True)
	print("Prediction:",clf.predict(predictQ))
	print("Actual answer",predictA,"\n")
	# Display the actual image
	displayImage(data,example)
	"""------------------------------Testing Function----------------------------"""
	# Actual testing on residual samples (=samples not used for training)
	print("Testing...")
	correct = 0
	wrong = 0
	j = example+1
	while j < getIndexMax(data):
	q = returnDataToUse(data,j,training=True,last=True)
	a = returnDataToUse(data,j,answers=True,last=True)
	p = clf.predict(q)
	if a == p[0]:
	correct += 1
	else:
	wrong += 1
	j += 1
	print("Statistics, correct answers:",correct/(correct+wrong))

view raw handwrittenDigits.py hosted with ❤ by GitHub

This time I got 89% success rate! Pretty good I guess! I wonder whether I could train Python to recognize other things, maybe faces or other! Well first of all I have to figure out how to convert a picture into readable numpy arrays. Readable for Python of course!! If you have any suggestion please do leave a comment!

	# >>> ================================ RESTART ================================
	# >>>
	# Maximum index: 1593
	#
	# Training...
	# Prediction: [9]
	# Actual answer 9
	#
	# [0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
	# Number to guess: 9
	#
	# Testing...
	# Statistics, correct answers: 0.8918918918918919
	# >>>

view raw handwrittenDigitsOut.py hosted with ❤ by GitHub

Here below is the citation of the source where I found the dataset “Semeion Handwritten Digits Data Set”:

Bache, K. & Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science.

and

Semeion Research Center of Sciences of Communication, via Sersale 117, 00128 Rome, Italy
Tattile Via Gaetano Donizetti, 1-3-5,25030 Mairano (Brescia), Italy.

 

Hope this was interesting!

Poker hands recognition (Machine Learning)

A few months ago I downloaded the scikit-learn package for Python, for those of you who might not be aware, scikit-learn is a powerful yet very simple package useful to apply machine learning. Basically they give you “all” the algorithms you may need and you “only” have to get the data and make it ready to feed into the scikit-learn algorithm.

Anyway, I only recently had time to check it out and write some code, furthermore only recently I found a great site full of sample datasets dedicated to machine learning (check it out here). Since the data from this great site is more or less in the right shape for being ready to import (.txt files or similar) the only task left to the programmer is to put it into numpy arrays.

This is one of my first attempt at generating a “real” machine learning program.

I won’t post the pre-settings script since it is pretty boring, instead I’ll briefly describe the pre-setting: basically, using a simple script I’ve splitted the training samples from the database in two .txt files trainingP.txt and target.txt respectively. TrainingP.txt contains the questions (the hands of poker) and target.txt contains the answers (=the score, that’s to say poker, full house etc… coded into numbers according to the description given in the database description file).

Below is the ml script: it is composed of 3 parts: setting, training and testing

setting: import the training sets and fit them into numpy arrays of the correct shape

training: fit the data into the model using scikit-learn

testing: test the algorithm and check what it has learned so far! Some statistics are printed, for reference. Check the results below!

	import numpy as np
	import sklearn
	import os
	# Import of support vector machine (svm)
	from sklearn import svm
	"""---------------------------- SETTING -----------------------------------"""
	# Changing current working directory
	os.chdir("C:\\")
	# Set up ML and scikit-learn
	clf = svm.SVC(gamma=0.01, C=100)
	# Training data (digits: questions, target: answers)
	digits = open("trainingP.txt").read().split("\n")
	target = open("target.txt").read().split("\n")
	# Initializing numpy arrays
	digitsD = np.array(np.fromstring(digits[0],sep=","))
	# Completing numpy array for digits
	for line in digits[1:]:
	try:
	sample = np.fromstring(line,sep=",")
	digitsD = np.vstack((digitsD,sample))
	except:
	pass
	# Converts strings into integers (not necessary if you
	# choose option 1 below)
	t = [int(i) for i in target]
	# 1 : All strings
	# 2 : All integers
	t = np.array(target) # 1
	#t = np.array(t) # 2
	"""---------------------------- TRAINING -----------------------------------"""
	# Model fitting (training)
	x,y = digitsD,t
	print("Training...")
	clf.fit(x,y)
	# Answer: 1
	print("Prediction:",clf.predict(np.array([3,11,1,9,2,5,3,5,4,1])))
	"""---------------------------- TESTING -----------------------------------"""
	# Load questions and answers
	rawQuestions = open("testQuestions.txt").read().split("\n")
	rawAnswers = open("Answers.txt").read().split("\n")
	# Initialize numpy arrays
	questions = np.array(np.fromstring(rawQuestions[0],sep=","))
	answers = np.array(rawAnswers)
	for q in rawQuestions[1:]:
	try:
	q1 = np.fromstring(q,sep=",")
	questions = np.vstack((questions,q1))
	except:
	pass
	# Testing function
	def tryAnswer(qu,an):
	machineAnswer = clf.predict(qu)
	if machineAnswer == an:
	#print("Correct!")
	return True
	else:
	#print("Wrong!")
	return False
	print("Testing...")
	# Start the test
	j = 0
	right = 0
	wrong = 0
	while j < 10000:
	result = tryAnswer(questions[j],answers[j])
	if result:
	right += 1
	else:
	wrong += 1
	j += 1
	# Print some statistics
	print("Statistics, correct rate:",right/(right+wrong))

view raw p_hands_recog.py hosted with ❤ by GitHub

So far the best score is 56.4%, I wonder if I did everything correctly! Anyway soon I will post a script with better results and achievements!

	# >>> ================================ RESTART ================================
	# >>>
	# Training...
	# Prediction: ['0']
	# Testing...
	# Statistics, correct rate: 0.564
	# >>>

view raw p_hands_recog2.py hosted with ❤ by GitHub

Below is the citation of the source of the database, according to their citation policy.

The name of the dataset is Poker Hand and it is from:

Bache, K. & Lichman, M. (2013). UCI Machine Learning Repository[http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science.

Fluid dynamics: pressure drop modelling

Likely the last post of the year, on a rather intriguing subject: fluid dynamics.

Recently I’ve been asked to model a simulation for a series of pipes through which water flows out of a reservoir. The task seems quite easy at first, however many factors must be taken into account. Pressure drop due to friction is one of the things that need to be addressed. Therefore, I documented myself a little on the subject and developed a simple model to take friction into account. Pressure drop depends both by the speed of the water and the pipes diameter.

Here is a short piece of code to calculate pressure drop. Note that we are using the average speed of the water. For each average speed we can draw a line that shows the pressure drop versus the pipe’s diameter.

Considering the fact that this was my first attempt I think the results are pretty good when compared with an actual diagram such as this (note that on this graph you’ll need to multiply pressure drop by 10 to compare pears to pears and apples to apples):

Here is the code I used. I should mention, I had a hard time finding out how to do this and furthermore I couldn’t find a comprehensive piece lecture where fluid dynamics is analysed in depth. Online there is a lot of material on Bernoulli’s equation and Torricelli’s special case but it is definitely harder to find out how to model real fluids, friction and fluid dynamics in general. Hope someone will upload some good material.

	from matplotlib import pyplot as plt
	import math
	# Density of the water
	rho = 1000
	# Kinematic viscosity of water
	nu = 1.01 * math.pow(10,-6)
	# Returns pressure drop Pa/m
	# Reynolds' number tells us which
	# expression to use in the two different scenarios
	# v velocity of the fluid
	# D pipes' diameter
	def pressureDrop(v,D):
	Re = (v*D)/nu
	if Re < 2000:
	print("Laminar flow",v,D)
	Fa = 64/Re
	else:
	print("Turbolent flow")
	Fa = 0.316 * math.pow(Re,-0.25)
	r = Fa*(1/D)*rho*math.pow(v,2)/2
	return r
	pressDropv = {}
	diaPipe = []
	veloc = []
	v = 0.5
	d = 0.05
	while v < 2.5:
	veloc.append(v)
	v += 0.5
	while d < 0.100:
	diaPipe.append(d)
	d += 0.01
	pressDropd = {}
	pressDropv = {}
	for i in veloc:
	pressDropv[i] = list()
	for j in diaPipe:
	pressDropv[i].append(pressureDrop(i,j))
	def plotAllD():
	try:
	x = diaPipe
	plt.title("Pressure drops given a fixed fluid velocity")
	plt.ylabel("Pressure drop Pa/m")
	plt.xlabel("Pipes' diameter m")
	for v in veloc:
	y = pressDropv[v]
	plt.plot(x,y,label="Fluid velocity "+str(v)+" m/s")
	plt.legend(loc="best", fontsize="small")
	plt.grid()
	plt.show()
	except Exception as e:
	print(e)
	plotAllD()

view raw pressure_drop.py hosted with ❤ by GitHub

And here comes the whole Solar System! ;)

Modelling Sun, Earth and Moon only? Naah that’s too naive. Let’s add some other planets!
The code on my previous article was rather ugly. It did not follow an OOP style and it was mostly non re-usable. Aside from the content of the piece of code I posted in this article, I’m quite satisfied with how I set it: a class for planets and a class to handle time. Sure, perhaps this is not the smartest way of using a class, however for someone like me, who programs mainly as a hobby, that’s fine I guess.
This piece of code models and simulates the whole Solar System (without the Sun, you can add it though!) and gives you, aside from the code, a hint on how slow/fast planets are orbiting around the Sun. note that I use large intervals of time (10000 seconds or similar). This model for sure can be improved. For instance the colours could be improved and, jokes aside, perhaps elliptical orbits could be implemented instead of circular. To do that I still need to work on my knowledge of gravity and its implementation in Python.


Anyway I’ll keep you posted. That’s it for now, here’s the code of the model and the YouTube video I made (but couldn’t load in the post, not sure why)

[EDIT]: I've finally made it and load the video here too! Enjoy

	import math
	from bigfloat import *
	import matplotlib.pyplot as plt
	from visual import *
	# A class to handle the time ranges
	class timeHoursSeconds(object):
	def __init__(self,s,h,d,y):
	self.s = s
	self.h = h
	self.d = d
	self.y = y
	def fromStoHours(self):
	h = self.s/60/60
	return h
	def fromStoDays(self):
	d = self.s/60/60/24
	return d
	def fromStoYears(self):
	y = self.s/60/60/24/365
	return y
	def fromDaysToS(self):
	s = self.d*24*60*60
	return s
	def fromDaysToH(self):
	h = self.d * 24
	return h
	def fromDaysToY(self):
	y = self.d/365
	return y
	class planet(object):
	G = 6.67 * math.pow(10,-11)
	sunM = 1.989 * math.pow(10,30)
	eaM = 5.973 * math.pow(10,24)
	RTL = 384400000
	def __init__(self,name,mass,RS,theta0,radius):
	self.name = name
	self.mass = mass
	self.RS = RS
	self.theta0 = theta0
	self.radius = radius
	def gravitationalForce(self,m2=1):
	if m2 ==1:
	f = self.G * (self.mass*self.sunM)/math.pow(self.RS,2)
	else:
	f = self.G * (self.mass*self.eaM)/math.pow(self.RTL,2)
	return f
	def angularVelocity(self,m2=1):
	w = math.sqrt(self.gravitationalForce(m2=m2)/(self.mass*self.RS))
	return w
	def velocity(self,m2=1):
	v = self.angularVelocity(m2=1) * self.RS
	return v
	def angularPosition(self,t,m2=1):
	theta = self.theta0 + self.angularVelocity(m2=m2) * t
	return theta
	def varAngularPosition(self,t,dt,m2=1):
	dtheta = self.angularPosition(t+dt,m2=m2)-self.angularPosition(t,m2=m2)
	return dtheta
	def periodAroundSun(self,m2=1):
	p = timeHoursSeconds(2*math.pi/self.angularVelocity(m2=m2),0,0,0)
	return p
	def picture(self,x,y,z,col,trail):
	if col == 1:
	return sphere(pos=vector(x,y,z),color=color.red,radius=self.radius,make_trail=trail)
	elif col == 2:
	return sphere(pos=vector(x,y,z),color=color.blue,radius=self.radius,make_trail=trail)
	elif col == 3:
	return sphere(pos=vector(x,y,z),color=color.green,radius=self.radius,make_trail=trail)
	elif col == 4:
	return sphere(pos=vector(x,y,z),color=color.cyan,radius=self.radius,make_trail=trail)
	elif col == 5:
	return sphere(pos=vector(x,y,z),color=color.yellow,radius=self.radius,make_trail=trail)
	else:
	return sphere(pos=vector(x,y,z),color=color.white,radius=self.radius,make_trail=trail)
	mercury = planet("Mercury",3.302 * math.pow(10,23),57910000000,0,0.3)
	venus = planet("Venus",4.8685 * math.pow(10,24),108200000000,0,0.4)
	earth = planet("Earth",5.973 * math.pow(10,24),149600000000,0,0.5)
	# As for the Moon, input Earth-Moon distance
	moon = planet("Moon",7.347 * math.pow(10,22),384400000,0,0.2)
	mars = planet("Mars",6.4185 * math.pow(10,23),227900000000,0,0.45)
	jupiter = planet("Jupiter",1.8986 * math.pow(10,27),778500000000,0,.8)
	saturn = planet("Saturn",5.6846 * math.pow(10,26),1433000000000,0,0.7)
	uranus = planet("Uranus",8.6832 * math.pow(10,25),2877000000000,0,0.6)
	neptune = planet("Neptune",1.0243 * math.pow(10,26),4503000000000,0,0.6)
	# Simulation data
	years = timeHoursSeconds(0,0,3655,0)
	seconds = years.fromDaysToS()
	print("Years: ",years.y)
	print("Days: ",years.d)
	print("Seconds: ",seconds)
	t = 0
	dt = timeHoursSeconds(10000,0,0,0)
	# Planets
	merc = mercury.picture(1.5,0,0,1,True)
	ven = venus.picture(3,0,0,3,True)
	ea = earth.picture(5,0,0,2,True)
	mar = mars.picture(7,0,0,3,True)
	jup = jupiter.picture(9,0,0,5,True)
	sat = saturn.picture(11,0,0,6,True)
	ur = uranus.picture(13,0,0,3,True)
	nep = neptune.picture(15,0,0,2,True)
	planetsf = [merc,ven,ea,mar,jup,sat,ur,nep]
	planets = [mercury,venus,earth,mars,jupiter,saturn,uranus,neptune]
	# The Moon
	v = vector(0.9,0,0)
	mo = moon.picture(5+0.9,0,0,10,True)
	for k in planets:
	revp = k.periodAroundSun()
	print("Planet name: ",k.name)
	print(k.name," mass: ",k.mass," kg")
	print(k.name," distance from the sun: ",k.RS/1000," Km")
	print(k.name," angular velocity: ",k.angularVelocity()," rad/s")
	print(k.name," period around the sun: ",revp.fromStoYears()," terrestrial year/s")
	print("\n")
	# Our program
	while t < seconds:
	rate(50)
	for plan in range(len(planets)):
	planetsf[plan].pos = rotate(planetsf[plan].pos,angle=planets[plan].varAngularPosition(t,dt.s),axis=(0,0,1))
	v = rotate(v,angle=moon.varAngularPosition(t,dt.s,m2=2),axis=(0,0,1))
	mo.pos = ea.pos + v
	t += dt.s

view raw solarSystem.py hosted with ❤ by GitHub

Earth Moon system orbiting around the Sun and VPython

Hello everyone! A few days ago my friend, as a joke, bet me to reproduce the Moon orbiting around the Earth. Well, I was a little afraid the model could be a little slow using matplotlib and its animation functions. Fortunately I found out VPython, a great 3D package: click here to get to their homepage. VPython, in my opinion, provides one of the simplest yet most brilliant package for 3D graphics. It’s shockingly easy to learn: it took less than an hour for me to learn the basics and start doing something interesting. Furthermore if you’re used to vectors and working with them, VPython provides simple functions for vectors calculations.

Check the model below (Earth and Moon orbiting around the Sun). I also made a video:

	import math
	from visual import *
	# Data in units according to the International System of Units
	G = 6.67 * math.pow(10,-11)
	# Mass of the Earth
	ME = 5.973 * math.pow(10,24)
	# Mass of the Moon
	MM = 7.347 * math.pow(10,22)
	# Mass of the Sun
	MS = 1.989 * math.pow(10,30)
	# Radius Earth-Moon
	REM = 384400000
	# Radius Sun-Earth
	RSE = 149600000000
	# Force Earth-Moon
	FEM = G*(ME*MM)/math.pow(REM,2)
	# Force Earth-Sun
	FES = G*(MS*ME)/math.pow(RSE,2)
	# Angular velocity of the Moon with respect to the Earth (rad/s)
	wM = math.sqrt(FEM/(MM * REM))
	# Velocity v of the Moon (m/s)
	vM = wM*REM
	print("Angular velocity of the Moon with respect to the Earth: ",wM," rad/s")
	print("Velocity v of the Moon: ",vM/1000," km/s")
	# Angular velocity of the Earth with respect to the Sun(rad/s)
	wE = math.sqrt(FES/(ME * RSE))
	# Velocity v of the Earth (m/s)
	vE = wE*RSE
	print("Angular velocity of the Earth with respect to the Sun: ",wE," rad/s")
	print("Velocity v of the Earth: ",vE/1000," km/s")
	# Initial angular position
	theta0 = 0
	# Position at each time
	def positionMoon(t):
	theta = theta0 + wM * t
	return theta
	def positionEarth(t):
	theta = theta0 + wE * t
	return theta
	def fromDaysToS(d):
	s = d*24*60*60
	return s
	def fromStoDays(s):
	d = s/60/60/24
	return d
	def fromDaysToh(d):
	h = d * 24
	return h
	# Graphical parameters
	print("\nSimulation Earth-Moon-Sun motion\n")
	days = 365
	seconds = fromDaysToS(days)
	print("Days: ",days)
	print("Seconds: ",seconds)
	v = vector(0.5,0,0)
	E = sphere(pos=vector(3,0,0),color=color.cyan,radius=.3,make_trail=True)
	M = sphere(pos=E.pos+v,color=color.white,radius=0.1,make_trail=True)
	S = sphere(pos=vector(0,0,0),color=color.yellow,radius=1)
	t = 0
	thetaTerra1 = 0
	dt = 5000
	dthetaE = positionEarth(t+dt)- positionEarth(t)
	dthetaM = positionMoon(t+dt) - positionMoon(t)
	print("delta t:",dt,"seconds. Days:",fromStoDays(dt),"hours:",fromDaysToh(fromStoDays(dt)),sep=" ")
	print("Variation angular position of the Earth:",dthetaE,"rad/s that's to say",degrees(dthetaE),"degrees",sep=" ")
	print("Variation angular position of the Moon:",dthetaM,"rad/s that's to say",degrees(dthetaM),"degrees",sep=" ")
	while t < seconds:
	rate(50)
	thetaEarth = positionEarth(t+dt)- positionEarth(t)
	thetaMoon = positionMoon(t+dt) - positionMoon(t)
	# Rotation only around z axis (0,0,1)
	E.pos = rotate(E.pos,angle=thetaEarth,axis=(0,0,1))
	v = rotate(v,angle=thetaMoon,axis=(0,0,1))
	M.pos = E.pos + v
	t += dt

view raw earthMoon.py hosted with ❤ by GitHub

Hope this was entertaining! Next comes the whole Solar System!

Gini coefficient, concentration measurement: an implementation in R

Another subject we took in the statistics class was the Gini index.

Gini index or ratio or coefficient is used to calculate how much a certain transferable phenomenon such as income or stocks for instance, is concentrated.

For example, say you are evaluating a company and you’d like to know more about how the shares are divided among the shareholders. You could use Gini index for that!

I’ve calculated the index using R and random data you can download here. In case you’d like to know more about Gini index check here.

Here my simple R implementation of the index.

	# Load data
	tb <- read.table("C:\\b.txt",header=TRUE,sep=",")
	# Add 5 new columns for analysis purposes
	for(i in 1:5)
	{
	cbind(tb,0)
	}
	# Storing the number of rows and columns
	nRow <- nrow(tb)
	nCol <- ncol(tb)
	# Cumulative frequencies
	i <- 1
	totalF = sum(tb[,2])
	while(i <= nRow)
	{
	if(i==1)
	{
	tb[1,3] <- tb[1,2]
	tb[1,4] <- tb[1,2]/1000
	}else{
	tb[i,3] <- tb[i-1,3]+tb[i,2]
	tb[i,4] <- tb[i-1,3]/totalF + tb[i,2]/1000
	}
	i <- i + 1
	}
	i <- 1
	while(i<=nRow)
	{
	tb[i,5] <- tb[i,1]*tb[i,2]
	if(i==1)
	{
	tb[i,6] <- tb[i,5]
	}else{
	tb[i,6] <- tb[i-1,6]+tb[i,5]
	}
	i <- i + 1
	}
	i <- 1
	while(i <= nRow)
	{
	tb[i,7] <- tb[i,6]/sum(tb[,5])
	i = i +1
	}
	# Show and plot the data
	tb
	a <- c(0,1)
	b <- c(0,1)
	c <- c(0,tb[,4])
	d <- c(0,tb[,7])
	plot(a,b,main="Concentration",type="l",col="green",lwd=2)
	lines(c,d,type="b",col="red",ylab="Relative freq",xlab="Relative freq",lwd=2)
	# Calculate Gini's R concentration index
	getR <- function(mat)
	{
	R <- 0.5
	area = 0.5*tb[1,4]*tb[1,7]
	i <- 2
	while(i <= nRow)
	{
	area = area + 0.5*(tb[i,4]-tb[i-1,4])*(tb[i-1,7]+tb[i,7])
	i = i + 1
	}
	acmax <- (sum(tb[,2])-1)/(2*sum(tb[,2]))
	R <- (R - area)/acmax
	return(R)
	}
	# Print data
	paste("Concentration index R is: ",getR(tb)*100,"%")

view raw concMeas.R hosted with ❤ by GitHub

Here below are the results

It looks like the data I used shows a 24% concentration. Cool!

Pearson’s chi-squared test: a simple implementation in R (test of independence)

Hi everyone! Today I found my old statistics workbooks and start wondering what I could get out of them.

Statistics can look pretty boring when using only pen and paper, since many times you’re just making a lot of repetitive calculations. However, the results of those calculations might of course be interesting.

Person’s chi-squared test is a simple test, as my professor put it, one of the first tests you should be performing when analysing a double entry table. You might be asking yourself why. Well, the answer is that this test looks for connection between the two variables in the table. As you might know connection is different from dependence. Dependence is of course a stronger bond and kind of a “one way bond” whilst connection is sort of a “double way bond”. If there’s no or little connection, then you might want to change variables in play since there is nothing but little interest in performing further tests on the same set of data.

Check on Wikipedia for more information on the theory.

Here is the R code I used to implement the test on the raw data at the bottom of the page.

	# Load data
	tb <- read.table("C:\\a.txt",header=TRUE,sep=",",row.names=1)
	# Add a column and a row
	tb <- cbind(tb,0)
	tb <- rbind(tb,0)
	# Store data on the table
	nRow <- nrow(tb)
	nCol <- ncol(tb)
	# Add row totals
	rownames(tb)[nRow] <- "Row sum"
	i <- 1
	while(i <= nRow)
	{
	partialRowsum <- sum(tb[i,])
	tb[i,nCol] <- partialRowsum
	i <- i + 1
	}
	# Add column totals
	colnames(tb)[nCol] <- "Columns sum"
	i <- 1
	while(i <= nCol)
	{
	partialColsum <- sum(tb[i])
	tb[nRow,i] <- partialColsum
	i <- i + 1
	}
	# This function returns a stocastic matrix
	# conditioned on the given one.
	stoc <- function(mat)
	{
	stocIndipTable <- matrix(0,nRow,nCol)
	k<-1
	while(k<=nRow)
	{
	j <- 1
	while(j<=nCol)
	{
	if((j==nCol) || (k==nRow))
	{
	stocIndipTable[k,j] <- tb[k,j]
	j = j + 1
	}else{
	stocIndipTable[k,j] <- tb[k,nCol]*tb[nRow,j]/tb[nRow,nCol]
	j = j + 1
	}
	}
	k = k + 1
	}
	return(stocIndipTable)
	}
	# Contingency table
	ctgTable <- tb-stoc(tb)
	# This function calculates the normalized
	# Pearson's chi squared (not constrained)
	calculatePearson <- function(cgt)
	{
	rawP <- cgt^2/stoc(tb)
	pearson <- sum(rawP)
	maxP <- min((nRow-2),(nCol-2)) * tb[nRow,nCol]
	return(pearson/maxP)
	}
	# Print information on the table
	paste("This is a ",dim(tb)[1]," by ",dim(tb)[2]," matrix.")
	print(colnames(tb))
	print(rownames(tb))
	# Pearson's Chi squared
	p <- calculatePearson(ctgTable) * 100
	if(p>0)
	{
	print("There is connection between the two phenomena in the table")
	print(paste("Chi squared connection index: ",p," %"))
	}else{
	print("Check process! Chi squared negative!")
	}

view raw chisqrd.R hosted with ❤ by GitHub

Below is the output, it seems there’s a feeble connection between the two phenomena we studied.

The data I used for the simulation are available to be downloaded here.

Here again Java

Java, Java and Java again! Here is my last script in Java, which I wrote today.

The script is quite simple, it is composed of three classes and lets you handle some basic calculations with lines and parabolas on the real plane and calculate the area under the curve (definite integrals).

Optionally you can plot the data.

Here is the main class with a simple example:

	public class Principal
	{
	public static void main(String[] args)
	{
	// Let's generate our line and parabola
	Integral f1 = new Integral(1,2,3);
	System.out.println(f1.description());
	System.out.println("\nIntersections between line and parabola");
	f1.findIntersections();
	System.out.println("\nIntersections of line with x axis");
	f1.findIntersectionsWithXaxis("line");
	System.out.println("\nIntersections of parabola with x axis");
	f1.findIntersectionsWithXaxis("parabola");
	// Integral of the parabola
	System.out.println("\nDefinite integral of the parabola from 0 to 10:");
	System.out.println(f1.approxIntegral(0,10,100000,"parabola"));
	// Plot of the function
	f1.plotFunction("both",-10,10,20);
	}
	}

view raw first_program1.java hosted with ❤ by GitHub

And the integral class which handles all the background operations

	import javax.swing.JFrame;
	import org.math.plot.*;
	public class Integral
	{
	private double a;
	private double b;
	private double c;
	// Constructor of the class
	public Integral(double a, double b, double c)
	{
	this.a = a;
	this.b = b;
	this.c = c;
	}
	// Returs a string with a simple description of what this class can do
	public String description()
	{
	String content = "This class lets you generate lines and parabolas"
	+ "\nFurthermore you can:"
	+ "\nCalculate definite integrals"
	+ "\nEvaluate functions"
	+ "\nPlot functions"
	+ "\nOutput data (txt format)";
	return content;
	}
	// Lets you evaluate the line at x and optionally print the result (print=true)
	public double line(double x,boolean print)
	{
	double value = a*x + b;
	if(print)
	{
	System.out.println("Line: the value of y at x="+x+" is equal to: "+value);
	}
	return value;
	}
	// Same as above but for the parabola
	public double parabola(double x, boolean print)
	{
	double value = a*Math.pow(x,2) + b*x + c;
	if(print)
	{
	System.out.println("Parabola: the value of y at x="+x+" is equal to: "+value);
	}
	return value;
	}
	// Evaluate definite integral with a finite number of steps
	public double approxIntegral(double initialp, double finalp, double steps, String function)
	{
	double sum = 0;
	double interval = (finalp-initialp)/steps;
	if(function.equals("line"))
	{
	while(initialp < finalp)
	{
	sum = sum + interval*this.line(initialp,false);
	initialp += interval;
	}
	return sum;
	}else if(function.equals("parabola"))
	{
	while(initialp < finalp)
	{
	sum = sum + interval*this.parabola(initialp,false);
	initialp += interval;
	}
	return sum;
	}else
	{
	String error = "Error, please insert either of the following:\nparabola\nline";
	System.out.println(error);
	return 0;
	}
	}
	// Plot the functions
	public void plotFunction(String function,double from, double to,int npoints)
	{
	if((from >= to) || (from==to) || (npoints==0))
	{
	System.out.println("Check input please! Something is wrong.");
	return;
	}
	double step = Math.abs((to-from)/npoints);
	Plot2DPanel plot = new Plot2DPanel();
	double x[] = new double[npoints];
	double y[] = new double[npoints];
	if(function.equals("line"))
	{
	for(int i=0; i < npoints; i++ )
	{
	x[i] = from;
	y[i] = this.line(from,false);
	from += step;
	}
	plot.addLinePlot("Line 1",x,y);
	}else if(function.equals("parabola"))
	{
	for(int i=0; i < npoints; i++ )
	{
	x[i] = from;
	y[i] = this.parabola(from,false);
	from += step;
	}
	plot.addLinePlot("Parabola 1",x,y);
	}else if(function.equals("both"))
	{
	double y1[] = new double[npoints];
	for(int i=0; i < npoints; i++)
	{
	x[i] = from;
	y[i] = this.parabola(from, false);
	y1[i] = this.line(from, false);
	from += step;
	}
	plot.addLinePlot("Parabola 1",x,y);
	plot.addLinePlot("Line1",x,y1);
	}else
	{
	System.out.println("Please check your input, function's name seems not ok.");
	return;
	}
	// Plot name and frame settings
	plot.setAxisLabel(0,"X");
	plot.setAxisLabel(1,"Y");
	JFrame frame = new JFrame("Graph");
	frame.setContentPane(plot);
	frame.setVisible(true);
	return;
	}
	// Find intersections between line and parabola
	void findIntersections()
	{
	double delta = Math.pow(this.b-this.a, 2)-4*this.a*(this.c-this.b);
	if(delta < 0)
	{
	System.out.println("Delta is negative, no solutions!");
	}else
	{
	double x1 = (this.a - this.b + Math.sqrt(delta))/(2*this.a);
	double x2 = (this.a - this.b - Math.sqrt(delta))/(2*this.a);
	System.out.println("Intersections at:"
	+ "\nx1: "+x1 + "y1: "+this.parabola(x1,false)
	+ "\nx2: "+x2 + "y2: "+this.parabola(x2,false));
	}
	}
	// Self explanatory ;)
	void findIntersectionsWithXaxis(String function)
	{
	if(function.equals("parabola"))
	{
	if(this.c==0)
	{
	double x2 = -this.b/this.a;
	System.out.println("x1: 0 y1: 0");
	System.out.println("x2: "+x2+" y2: "+this.parabola(x2,false));
	}
	else
	{
	double delta = Math.pow(this.b, 2)-4*this.a*this.c;
	if(delta < 0)
	{
	System.out.println("Negative delta, no solutions!");
	return;
	}else
	{
	double x1 = (-this.b + Math.sqrt(delta))/(2*this.a);
	double x2 = (-this.b - Math.sqrt(delta))/(2*this.a);
	System.out.println("x1: "+x2+"y1: "+this.parabola(x1,false));
	System.out.println("x2: "+x2+"y2: "+this.parabola(x2,false));
	}
	}
	}else if(function.equals("line"))
	{
	double x1 = -this.b/this.a;
	System.out.println("x1: "+x1+"y1: "+this.line(x1,false));
	}
	}
	}

view raw first_program2.java hosted with ❤ by GitHub

Optionally, should you want to output something to a .txt file, here is the class that can help in doing that:

	import java.io.BufferedWriter;
	import java.io.File;
	import java.io.FileWriter;
	import java.io.IOException;
	public class GetOutData
	{
	private String content;
	// if you need to print a string
	public GetOutData(String output)
	{
	this.content = output;
	}
	public void writeString(String path)
	{
	try
	{
	File file = new File(path);
	if(!file.exists())
	{
	file.createNewFile();
	}
	FileWriter fw = new FileWriter(file.getAbsoluteFile());
	BufferedWriter bw = new BufferedWriter(fw);
	bw.write(this.content);
	bw.close();
	}catch(IOException e)
	{
	e.printStackTrace();
	}
	}
	}

view raw first_program3.java hosted with ❤ by GitHub

Bye until the next project!

Java again, in English this time though! ;)

As i said earlier in my article in German, I studied Java only once in a while, since I had so many other things to do and study. Anyway here is another bit of code I’ve written in order to practice what I’ve learnt.

This script is a simple program to keep track of your inventory. Frankly speaking, an excel spreadsheet would have been faster and easier but where’s the fun in it?

There are three classes: the main class, the products class and the inventory class. In the products class lies all the code related to the object “product” while in the inventory class you can find all the code needed to run the inventory. Finally the main class executes the main program and (as it should be I guess).

For the sole purpose of exercising myself with Java, I tried to generate a inventory of food with three items. However you can easily add as many inventories as you like and figure out a way to speed up the process of adding items to each inventory.

Perhaps in the future I’ll add the option to print out a bar chart or something similar as soon as I have time to study Java. As for now, I find it quite interesting and demanding since it asks every time for the type of the data you are going to work with and the arrays are different from python’s lists. That can be a small issue for someone who got used to Python as myself. Furthermore Java libraries are hugely vast and that can be overwhelming at the beginning. Nonetheless I hope to get better at it soon!

Here is the products class

public class Products 
{
	private String name;
	private int id;
	private int quantity;
	private double price;
	
	//Constructor of the class
	public Products(String name,int id, int quantity, double d)
	{
		this.name = name;
		this.id = id;
		this.quantity = quantity;
		this.price = d;
	}
	
	//Quantity setter and update method
	public void changeQuantity(int q,String s)
	{
		if(s.equals("add"))
		{
			this.quantity += q;
		}else if(s.equals("subtract"))
		{
			this.quantity -= q;
		}
		System.out.println("Warehouse stock updated succesfully!");
	}
	
	//Price setter
	public void changePrice(double p)
	{
		this.price = p;
		System.out.println("Price changed succesfully!");
	}
	
	//Get all the info on the product
	public void getInfoProduct()
	{
		System.out.println("As for product "+this.name+" with id: "+this.id);
		System.out.println("Quantity available: "+this.quantity);
		System.out.println("Price "+this.price);
		System.out.println("Total value in stock: " + this.getTotalValue());
	}
	
	//This function returns the total value of the stock for the product
	public double getTotalValue()
	{
		double totalValue = this.price*this.quantity;
		return totalValue;
	}
	
}

The main class

public class mainClass 
{

	public static void main(String[] args) 
	{	
		//we create the products
		Products bread = new Products("Bread",0,10,0.5);
		Products ooil = new Products("Olive oil",1,20,4.00);
		Products oranges = new Products("Oranges",2,5,2.50);
		
		//we create the inventory and add products to it
		Inventory inventory1 = new Inventory("Inventory 1");
		inventory1.addProduct(bread);
		inventory1.addProduct(ooil);
		inventory1.addProduct(oranges);
		
		//We print info on the products
		bread.getInfoProduct();
		ooil.getInfoProduct();
		oranges.getInfoProduct();
		
		//And the total value of the inventory
		System.out.println("\nTotal value of the inventory is "+inventory1.getInventoryValue());
	}
}

Finally, the inventory class

//Dynamics arrays. I found them similar to those in Python
import java.util.ArrayList;
import java.util.List;

public class Inventory 
{
	private String name;
	//Products[] productsInStock = new Products[]{};
	List<Products> productsInStock = new ArrayList<Products>();
	private int totalItems = 0;
	
	//Constructor of the class
	public Inventory(String name)
	{
		this.name = name;
	}	
	
	//Add a product object to the inventory array
	public void addProduct(Products p)
	{
		productsInStock.add(p);
		this.totalItems += 1;
		System.out.println("Product added to inventory "+this.name+" succesfully!");
	}
	
	//Get the total value of the inventory
	public double getInventoryValue()
	{
		double valueToReturn=0;
		int lenAr = this.totalItems;
		for(int i=0;i<lenAr;i++)
		{
			valueToReturn += this.productsInStock.get(i).getTotalValue();
		}
		return valueToReturn;
	}
}

Here below is the output

Product added to inventory Inventory 1 succesfully!
Product added to inventory Inventory 1 succesfully!
Product added to inventory Inventory 1 succesfully!
As for product Bread with id: 0
Quantity available: 10
Price 0.5
Total value in stock: 5.0
As for product Olive oil with id: 1
Quantity available: 20
Price 4.0
Total value in stock: 80.0
As for product Oranges with id: 2
Quantity available: 5
Price 2.5
Total value in stock: 12.5

Total value of the inventory is 97.5

Hope this was interesting.

Es ist endlich Zeit mit Java zu spielen

Heute hatte ich Lust ein bisschen zum Thema Java zu schreiben, und der Artikel auf Deutsch zu schreiben um meine Kenntnisse der deutschen Sprache zu verbessern.

Erst muss ich mich entschuldigen, weil es im Artikel Fehler seien könnten.
Ich studiere Deutsch seit zwei Jahren als Hobby, aber ich weiss noch nicht, wenn mein Deutsch gut genug ist, um einen kleinen Artikel zu Schreiben. Auf jeden Fall, wollte ich es ausprobieren, weil man sich nur mit Praxis verbessert.

Obwohl man in einem kleinen Artikel wie meiner nicht so viele Wörter braucht, kann es schwierig sein gut Deutsch zu schreiben. Aber ich denke, dass die komplizierte Grammatik der deutsche Sprache ein sehr attraktiver Aspekt sein kann.

Ich habe genug über die Sprache geredet, jetzt reden wir über Java.

Ich bin kein Experte von Java, ich habe seit dem letzten Sommer Java studiert aber nicht jeden Tag sondern nur ab und zu, wenn Zeit da war. Jedenfalls ist hier eines der ersten Projekte. Es ist ein sehr kurzes Skript, das euch eine Grafik zeigen kann.

Um das Skript zu machen, habe ich Eclipse benutzt. Auf dieser website kann man Eclipse finden.

Hier ist das Skript

	import javax.swing.JFrame;
	import org.math.plot.*;
	import java.util.ArrayList;
	import java.util.List;
	public class Principale
	{
	public static void main(String[] args)
	{
	UniformlyAcceleratedMotion a = new UniformlyAcceleratedMotion(0,10,2);
	System.out.println(a.getPosition(2));
	System.out.println(a.getSpeed(2));
	a.printGeneralSettings();
	// List position = new ArrayList();
	double space[];
	space = new double[10];
	double speedd[];
	speedd = new double[10];
	double time[];
	time = new double[10];
	for(int i=0;i<10;i++)
	{
	space[i] = a.getPosition(i);
	speedd[i] = a.getSpeed(i);
	time[i]=i;
	}
	// plot2DPanel object
	Plot2DPanel plot = new Plot2DPanel();
	// add line
	plot.addLinePlot("Space",time,space);
	plot.addLinePlot("and Speed",time,speedd);
	// Plot name and frame settings
	JFrame frame = new JFrame("Uniformly accelerated motion");
	frame.setContentPane(plot);
	frame.setVisible(true);
	}
	}

view raw java_example1.java hosted with ❤ by GitHub

Das Skript benutzt das Class “UniformlyAcceleratedMotion”, das hier unter ist:

	public class UniformlyAcceleratedMotion
	{
	private double v0;
	private double x0;
	private double a;
	public UniformlyAcceleratedMotion(double x0, double v0, double a)
	{
	this.x0 = x0;
	this.v0 = v0;
	this.a = a;
	}
	public double getPosition(double t)
	{
	if(t <= 0)
	{
	return this.x0;
	}else
	{
	double position = this.x0 + this.v0 * t + 0.5*this.a * Math.pow(t,2);
	return position;
	}
	}
	public double getSpeed(double t)
	{
	if(t<=0)
	{
	return this.v0;
	}
	else
	{
	double speed = this.v0 + this.a * t;
	return speed;
	}
	}
	public void printGeneralSettings()
	{
	System.out.println("Initial position: "+this.x0 + " m");
	System.out.println("Initial speed: "+this.v0 + " m/s");
	System.out.println("Acceleration (constant): "+this.a + " m/s^2");
	}
	}

view raw java_example2.java hosted with ❤ by GitHub

Hier ist die Graphik:

Ich hoffe, dass dieser Artikel lustig war. Auf wiedersehen, bis zum nächsten Mal.
Vielen Dank!

A floating ball dropped in a water current

While writing the previous script on simple harmonic motion I found a funny simulation.
Here it is, do you remember what happens when you drop a ball in flowing water? (Actually the fluid simulated looks something more dense than water).

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	import random
	fig = plt.figure()
	ax1 = fig.add_subplot(1,1,1)
	def fun(x):
	if x!=0:
	return 2*np.sin(4/3*x+2)/x
	else:
	return 1
	counter = 1
	def animate(i):
	global x,y,counter
	counter += 0.2
	ax1.clear()
	plt.xlim(0,50)
	plt.ylim(-2,2)
	plt.plot([x for x in range(50)],[-0 for x in range(50)],lw=2,color="cyan")
	plt.plot(counter,fun(counter),color="red",marker="o",ms=20)
	ani = animation.FuncAnimation(fig,animate,interval=1)
	plt.show()

view raw ball_water.py hosted with ❤ by GitHub

Basic physics and Python: simple harmonic motion

Here is simple harmonic motion simulation with a spring and a bouncing ball.
Springs are a classic example of harmonic motion, on Wikipedia you can get a grasp of the basics. Among other assumption, in my simulation I’ve assumed an ideal spring and that there is no friction (and therefore the motion will not stop by itself) however, if you like, you can implement friction easily.

Here is the spring simulation

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	import random
	fig = plt.figure(figsize=(7,7))
	ax1 = fig.add_subplot(1,1,1)
	ax2 = fig.add_subplot(2,1,1)
	ax3 = fig.add_subplot(3,1,1)
	fig.subplots_adjust(hspace=.45)
	# Spring data
	k = 100
	m = 20
	w = np.sqrt(k/m)
	phi = 2
	A = 2
	period = 2*np.pi/w
	frequency = 1/period
	print("Period:",period,"seconds",sep=" ")
	print("Frequency:",frequency,"Hz",sep=" ")
	print("k: ",k,"N/m",sep=" ")
	def fun(t):
	global w,phi,A
	return A*np.sin(w*t+phi)
	def vel(t):
	global w,phi,A
	return A*w*np.cos(w*t+phi)
	def acceleration(t):
	global w,phi,A
	return -A*w**2*np.sin(w*t+phi)
	def position(x):
	x1 = x-1
	x2 = x+1
	y1 = 1
	y2 = -1
	p1 = [x1,y2]
	p2 = [x2,y2]
	p3 = [x2,y1]
	p4 = [x1,y1]
	return [p1,p2,p3,p4]
	counter = 0
	xt = [0]
	yt = [0]
	vy = [0]
	acy = [0]
	def animate(i):
	global counter, xt,yt,A,vy,acy
	ax3.clear()
	plt.subplot(311)
	# configure X axes
	plt.xlim(-3.5,3.5)
	# configure Y axes
	plt.ylim(-2,2)
	# labels
	plt.xlabel("x")
	plt.ylabel("y")
	plt.title("Motion of an ideal spring")
	p1 = [position(fun(counter))[0][0],position(fun(counter))[0][1]]
	p2 = [position(fun(counter))[1][0],position(fun(counter))[1][1]]
	p3 = [position(fun(counter))[2][0],position(fun(counter))[2][1]]
	p4 = [position(fun(counter))[3][0],position(fun(counter))[3][1]]
	x = [p1[0],p2[0],p3[0],p4[0],p1[0]]
	y = [p1[1],p2[1],p3[1],p4[1],p1[1]]
	linex = [-4,p1[0]]
	liney = [0,0]
	plt.plot(x,y,lw=5,color="blue")
	plt.plot(linex,liney,color="red",ls=":",lw=5)
	plt.subplot(312)
	xt.append(counter)
	vy.append(vel(counter))
	plt.title("Velocity")
	plt.xlim(0,15)
	plt.ylim(-A*w-0.5,A*w+0.5)
	plt.plot(xt,vy,lw=1,color="green")
	plt.plot([0,15],[0,0],lw=0.5,color="black")
	plt.subplot(313)
	acy.append(acceleration(counter))
	plt.title("Acceleration")
	plt.xlim(0,15)
	plt.ylim(-A*w**2-0.5,A*w**2+0.5)
	plt.plot(xt,acy,lw=1,color="green")
	plt.plot([0,15],[0,0],lw=0.5,color="black")
	counter += 0.1
	ani = animation.FuncAnimation(fig,animate,interval=2)
	plt.show()

view raw harmonicMotion1.py hosted with ❤ by GitHub

Below is the bouncing ball:

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	import random
	fig = plt.figure()
	ax1 = fig.add_subplot(1,1,1)
	def fun(x):
	return 2*np.sin(2*x+2)
	counter = 1
	def animate(i):
	global x,y,counter
	counter += 0.2
	ax1.clear()
	plt.xlim(0,50)
	plt.ylim(-3,3)
	plt.plot([x for x in range(50)],[-2.15 for x in range(50)],lw=2,color="black")
	plt.plot(counter,fun(counter),marker="o",ms=20)
	ani = animation.FuncAnimation(fig,animate,interval=1)
	plt.show()

view raw harmonicMotion2.py hosted with ❤ by GitHub

A simpler example of simple harmonic motion with a spring (video):

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	import random
	fig = plt.figure()
	ax1 = fig.add_subplot(1,1,1)
	def fun(t):
	return 2*np.sin(4/3*t+2)
	x = [0]
	y = [fun(x[0])]
	def position(x):
	x1 = x-1
	x2 = x+1
	y1 = 1
	y2 = -1
	p1 = [x1,y2]
	p2 = [x2,y2]
	p3 = [x2,y1]
	p4 = [x1,y1]
	return [p1,p2,p3,p4]
	counter = 0
	def animate(i):
	global x,y,counter
	ax1.clear()
	# configure X axes
	plt.xlim(-5,5)
	plt.xticks([-4,-3,-2,-1,0,1,2,3,4])
	# configure Y axes
	plt.ylim(-5,5)
	# labels
	plt.xlabel("x")
	plt.ylabel("y")
	plt.title("Simple harmonic motion")
	# plt.gca().add_line(line)
	p1 = [position(fun(counter))[0][0],position(fun(counter))[0][1]]
	p2 = [position(fun(counter))[1][0],position(fun(counter))[1][1]]
	p3 = [position(fun(counter))[2][0],position(fun(counter))[2][1]]
	p4 = [position(fun(counter))[3][0],position(fun(counter))[3][1]]
	x = [p1[0],p2[0],p3[0],p4[0],p1[0]]
	y = [p1[1],p2[1],p3[1],p4[1],p1[1]]
	linex = [-4,p1[0]]
	liney = [0,0]
	plt.plot(x,y,lw=5,color="blue")
	plt.plot(linex,liney,color="red",ls=":",lw=5)
	counter += 0.2
	ani = animation.FuncAnimation(fig,animate,interval=1)
	plt.show()

view raw harmonicMotion3.py hosted with ❤ by GitHub

Hope this was interesting.

Animated graphs with matplotlib

Recently I have had not so much time to dedicate to my blog, however, today I had some spare time and decided to learn how to code animated graphs with Python and matplotlib.
Apparently it’s really simple, however a little bit of practice is needed, here below are three pieces of code where I coded and plotted the graph of (many) branch of parabolas, a random walk, and a simulated stock (a random walk again!).
Here is the half parabolas

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	fig = plt.figure()
	ax1 = fig.add_subplot(1,1,1)
	x = [i for i in range(20)]
	y = [np.sqrt(i**3) for i in x]
	def animate(i):
	global x,y
	r = np.random.normal(0,1,1)
	xar = x
	yar = y + r*y
	#ax1.clear()
	ax1.plot(xar,yar,marker='o')
	ani = animation.FuncAnimation(fig,animate,interval=1000)
	plt.show()

view raw animatedGraphs1.py hosted with ❤ by GitHub

Below is the random walk and the video:

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	fig = plt.figure()
	ax1 = fig.add_subplot(1,1,1)
	def randomWalk():
	r = np.random.normal(0,1,1)
	if r < 0:
	return -1
	elif r > 0:
	return 1
	else:
	return 0
	counter = 0
	start = [0]
	t = [0]
	def animate(i):
	global t,u,sd,counter
	x = t
	y = start
	counter += 1
	x.append(counter)
	y.append(start[counter-1] + randomWalk())
	plt.plot(x[counter-1],y[counter-1],marker='o',color="r")
	ani = animation.FuncAnimation(fig,animate,interval=1)
	plt.show()

view raw animatedGraphs2.py hosted with ❤ by GitHub

video of the random walk: https://www.youtube.com/watch?v=3pAHwt1ioz4

and finally the stock simulation:

	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	import random
	fig = plt.figure()
	ax1 = fig.add_subplot(1,1,1)
	def getNewPrice(s,mean,stdev):
	r = np.random.normal(0,1,1)
	priceToday = s + s*(mean/255+stdev/np.sqrt(225)*r)
	return priceToday
	z = np.random.normal(0,1,255)
	u = 0.1
	sd = 0.3
	counter = 0
	price = [100]
	t = [0]
	def animate(i):
	global t,u,sd,counter
	x = t
	y = price
	counter += 1
	x.append(counter)
	y.append(getNewPrice(price[counter-1],u,sd))
	ax1.clear()
	plt.plot(x,y,color="blue")
	ani = animation.FuncAnimation(fig,animate,interval=50)
	plt.show()

view raw animatedGraphs3.py hosted with ❤ by GitHub

Hope this was interesting .