How to fit data to a normal distribution using MLE and Python

MLE, distribution fittings and model calibrating are for sure fascinating topics. Furthermore, from the outside, they might appear to be rocket science. As far I'm concerned, when I did not know what MLE was and what you actually do when trying to fit data to a distribution, all these tecniques did looked exactly like rocket science.
They are not that much complicated though. MLE is a technique that enables you to estimate the parameters of a certain random variable given only a sample by generating a distribution which makes the observed results the most likely to have occurred. Distribution fittings, as far as I know, is the process of actually calibrating the parameters to fit the distribution to a series of observed data.

Let's see an example of MLE and distribution fittings with Python. You need to have installed scipy, numpy and matplotlib in order to perform this although I believe this is not the only way possible. For some reason that I ignore, the methods in scipy.stats related to the normal distribution use loc to indicate the mean and scale to indicate the standard deviation. I maybe can grasp why use "scale" to indicate the stdv however I really do not get "loc" I do not understand why... If you know that, please leave a comment.

The result should look somewhat like this:

Hope this was useful.

Estimating the area under a curve using random points in R

Even though this post is not the first one on this blog, the code I used to estimate the area under the curve x^2 and posted below is one of my first attempts to code something "complicated" with R. It has been more than 5 months since I wrote it, it does look a bit clunky now. Indeed. However I believe it is fun to keep an eye on our own written code and see how it has evolved. Even though the code is lengthy, the final result looks nice.

 

Here is the result you should get:

Simulate data with R

Last semester I was attending a boring class, even though the professor was really clever, he was always bouncing around the main theme and never got straight to the point. While thinking about everything but the class, I had an idea: when you are given a set of data, say X and Y, you can easily compute a linear regression model, e.g. the regression line, and find out information on the data. Now, you will also find information on the error that the linear model made in predicting the data. By finding out the distribution of the error you can somehow simulate data similar to the original, from the regression line, by simply adding a random error (whose distribution is known) to the predicted data.
Furthermore, we know from the regression line that the expected error is 0.

Here is the code to implement this idea in R. You can get the data to work on in the bottom of the page.


The result should look something like this: In blue the actual data and in red the simulated one.

Hope this was useful, if you know the name of this method, please leave a comment and let me know. Click here to get the data I used.

Copula functions in R

A copula function is an application which "couples" (joins) a multivariate distribution to its univariate margins (marginal distributions).
Copula functions can be really helpful in building multivariate distributions given the marginals. Here is a fast introduction to copulas.


A copula C can be defined as follows:

where I is the interval [0,1].


Archimedean copulas are a particular class of copulas which can be built using a function phi known as the copula generator, from the following relation:

In this post, we are going to see the main formulas for using a particular Archimedean copula in R: the Gumbel copula.
The gumbel copula is built using the generator function below

and has the following expression

The package, available in R, which has some useful functions on the gumbel copula is called gumbel, you need to install it and then call it like this

library(gumbel)

Once you called it, here are some basic functions
here is the density function

#plot the density
x <- seq(.01, .99, length = 50)
y <- x
z <- outer(x, y, dgumbel, alpha=2)
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightgreen", ltheta = 100,xlab = "x",ticktype = "detailed", ylab = "y", zlab = "Density of the Gumbel copula")

and the cumulative distribution function (CDF)

z <- outer(x, y, pgumbel, alpha=2)
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightgreen",ltheta = 100, ticktype = "detailed",xlab = "u", ylab = "v", zlab = "Cumulative distribution function")

Finally, we are going to take a look at the random number generating function. The range of dependece simulated by the Gumbel copula, depends only on the parameter theta. As theta increases, so does the dependence between observations. As you may have noticed, if theta is equal to 1 (theta is in [1,Inf) for the Gumbel copula), then we fall back in the independece case here below

#we simulate 2000 observations with theta = 1
r_matrix <- t(rgumbel(2000,1))plot(r_matrix[1,], r_matrix[2,], col="blue", main="Gumbel, independence case")

When we increase theta we obtain a different result, as expected

#we simulate 2000 observations with theta = 2
r_matrix <- (rgumbel(2000,3))plot(r_matrix[1,], r_matrix[2,], col="blue", main="Gumbel, Positive dependence")

you can see that Gumbel copula can be used to simulate positive and asymetric dependence, in fact the correlation seems to be higher on larger values. Below an example for theta = 3

On YouTube, I uploaded a simple animation created with R and windows movie maker you can watch it in the embedded video below

Hope this was useful.

Simulating random points with Python

Random, randomness and order are surely fascinating topics. In fact, even thought one might think random is a relatively easy idea, it is actually much more complicate than one might think.
This video is really interesting and if you are interested in randomness and random variables you should definetly check it out:https://www.youtube.com/watch?v=nAxEzxHkqyY


Computers can simulate random values if ordered to do so. Although they appear random these values are not really random. In fact, they are generated through a certain procedure and can be replicated. This feature of course can be useful if you run some kind of simulation with random numbers and then you want to replicate your simulation with the very same "random" numbers. Remember the line set.seed(a)? If you run this line in Python before you run your script, it will generate random number in a certain way according to the parameter (a) which should be a number (I believe positive only are allowed). Next time when you run the script again, if you want to get the same "random" numbers, you just need to add the line set.seed(a) and make sure the parameter a is the same you used the last time.

Simulating values can be very useful, say for instance that you know the distribution of a random variable for example you know that each number on a dices has the same probability of showing. Now you can simulate that random variable. Python and R have many interesting functions to generate random numbers. In future I am going to make a comparison post between Python and R. Today I am going to show you just some basic random generating functions in Python.

Random uniform distribution:
The probability density function of the uniform distribution is the following

The probability mass is uniformely distributed in the interval [a,b]. In Python the function random.uniform(a,b) generates a random value x, a < x < b. Let's see an example where a = 0 and b = 1:

import random

from matplotlib import pyplot as plt

import pylab



#random uniform sample

sample1 = []

sample2 = []



i = 0

while i < 10000:
    sample1.append(random.uniform(0,1))
    sample2.append(random.uniform(0,1))
    i += 1

plot = plt.plot(sample1,sample2, "bo")
    
plt.show()

Here is the result if we plot the two random vectors we have generated above:

In fact it could be fun to combine each random number generating function with the others and plot the result. Here some examples:

sample1 = []
sample2 = []

i = 0
while i < 2000:
    sample1.append(random.gauss(0,1))
    sample2.append(random.uniform(3,1))
    i += 1

plot = plt.plot(sample1,sample2, "bo")
plt.show()

and the result:

In this example we combine a uniform distribution with a gaussian one. Let's see a final example with gammavariate distribution and gaussian

sample1 = []
sample2 = []

i = 0
while i < 2000:
    sample1.append(random.gauss(0,1))
    sample2.append(random.gammavariate(2,1))
    i += 1

plot = plt.plot(sample1,sample2, "bo")
plt.show()

and the result:

With this method you can more or less simulate different behaviour of random variables by combining different distributions. Perhaps I'll post more examples when I'll write about the random module in Python and R.

Another tool which is useful to study random variables and their joint behaviour are copulas. A copula is a function which joints together many CDF and returns the joint CDF (Cumulative Distribution Function). It enables you to express the joint cumulative distribution of two or more random variables as a function of their marginals. As far as I know, only R has some functions for copulas, perhaps I'll make a post on it in the future.

You can check out this wikipedia page for more information on copulas: http://en.wikipedia.org/wiki/Copula_%28probability_theory%29

The Monty Hall problem

The Monty Hall problem is a famous game which was played in the television show "Let's make a deal".

The game goes like this:

There are three doors, behind each door there is either a goat or an amazing sportcar. The contestant wins if they guess where the car is. There are in total 2 goats and the car. Therefore the initial chance of choosing the car is 1/3.

The host asks the contestant to pick a door. Once a door has been chose, the host, who knows where the goats and the car are, opens a door behind which there is a goat and asks the contestant if they want to switch door.

Now the question is: should the contestant change door or should they stay? Is there any statistical reason which could justify either choice?

It might not be that immediate to understand, however the optimal strategy is to change door. In this case, the probability of winning the car increases from 1/3 to 2/3. You can check this by analysing the favourable scenarios over the possible scenarios. However, should we stick with this or should we make a simulation to test this statement out? Let's go for the simulation with R.

Here are the results:

The simulation confirms that, on the long run, the odds are more favourable if the contestant decides to switch door. Hope this was useful and interesting.
More on the Monty Hall problem: http://en.wikipedia.org/wiki/Monty_Hall_problem

A small sidenote on cx_Freeze, converting (GUI) scripts into exe and the annoying dos window

Hi everyone!

As you know, Python scripts need to be interpreted in order to run on your computer. Let's say that you want to give your brand new program to your friend who, for some reason, is reluctant to install Python and does not like using the command line or related things. Well then building a simple GUI (Graphical User Interface), with Tkinter for instance, is the right choice. However, we have not solved the user friendliness issue since even with a GUI, our script needs to be run with Python. The solution then, is using the module cx_Freeze to "convert" our script into an executable file (.exe) which can be run on windows.
I have been using cx_Freeze for just 3 months I guess, and I find it practical and useful. There are other modules which serve the same purpose, however some of them (py2exe for instance) are still not available for Python 3. More information on http://cx-freeze.sourceforge.net/

Here are the main steps to convert your .py program into an .exe:
1.Check that your script works fine
2.Create a setup.py file as below
3.Put setup.py and your script in the folder pythonxx (python33 if you have python 3.3)
4.Open the command line (cmd) and go to the pythonxx folder
5.Type in: python33\python.exe setup.py build
6.Press enter. Python will build your exe and put a "Build" folder in your python33 folder.
7.Check out the Build folder in python33 where the executable lies. There you go!


The setup file should contain this (sqrt.py is the script to converted):

from cx_Freeze import setup, Executable

setup(
    name = "sqrt",
    version = "0.1",
    description = "test",  #you can put here whatever you want
    executables = [Executable("sqrt.py")],
    )

Once you have had fun with your new executable script, if your script had a GUI interface, such as Tkinter, you will find an annoying thing: the dos window!!! It will not go away no matter what you do. It just comes up every time you run the .exe. The solution to this problem: you need to edit the setup.py file and reconvert the script. Here is the setup.py if your script has a GUI such as Tkinter:

import sys
from cx_Freeze import setup, Executable

base = None
if (sys.platform == "win32"):
    base = "Win32GUI"


setup(
    name = "sqrt",
    version = "0.1",
    description = "test",
    executables = [Executable("sqrt.py",base=base)],
    )

And your dos window will not appear every time you run the .exe. I had a hard time finding this and I hope it will be useful to someone.
Enjoy!

Physics with Python

Today I am going to built a class in Python to simulate a famous kind of motion: projectile motion!

I said a "Class" because while trying to find something to build a class on, I noticed that a projectile could fit perfectly. In fact, a projectile, or any other object which when launched, moves along a curved path under the action of gravity only, has some characteristics: speed, mass, etc...
Displaying all these characteristics through a class can be at the same time both fun and useful for practicing with classes.

The basic assumptions of this kind of motion are the following:
- There is no air friction
- The only significant force which acts on the projectile is gravity

Additional assumptions:
- Projectile starts at y = 0 (from ground)


If you'd like to gather more information on this kind of motion, check out wikipedia:
http://en.wikipedia.org/wiki/Projectile_motion



Here are the results:

Here is a video I made with the instance p and the last method:

It would be fun to implement air friction and then compare the data.

If you find incorrect physics terms please let me know or comment, I might not fully remember my physics classes!! :)

Approximating the value of pi with a Monte Carlo approach

Sometimes a value cannot be found so easily, this is the case for probabilities for example. When it comes to find the chances, say of a football team winning a certain match, you cannot use traditional probability approach, which is:

or, I should say, you could use it, but it can be misleading. For instance, the probability of winning for each team should be 50% even if the two teams are a Premier League team and a local team... is there anyone who believes this?
A better approach in these cases is simulating the entire process (i.e. the game) a big number of time. The bigger the number of simulations, the more precise the parameter(s) estimate. In case of a football match for example, you could just model all the players and their characteristics and simulate thousands of games. Even in this case you do not have a 100% confidence that your estimate is correct, since the initial parameters in the simulations, of course, affect the outcome and therefore a wrong initial parameter could lead to a wrong estimate. However, you can get interesting results using this approach, especially if you apply it to games of chance.

Another use of this method is the approximations of parameters such as pi, areas, and other values which are not easily computable in another way or that for some reason you do not want to compute in another way.

We can use R and Monte Carlo method to approximate the value of pi, and this is exactly what this post aims to.
By the way, approximating values using Monte Carlo method is one of my favourite techniques.

Perhaps in other posts I will write about Monte Carlo methods and games of chance.

The approach to calculate pi using Monte Carlo is pretty straightforward:
Say you have a circumference of radius 1, you can fit a square around it and easily calculate the area of the square.
Now, imagine to fill the square with thousands of dots. (You can do the same with a semi-circumference and it will be easier to code that, therefore we will opt for this solution).

The more dots you use, the more precise the approximation of pi.

Some of them will be inside the circumference and some of them will be outside of it. By calculating the proportion of inside dots and outside dots you can find the ratio of the circumference area to the square area. Now you can approximate the circumference area and reversing the formula for the area, obtain an approximation of pi.

In this example we use 100.000 points.


The output should be something like:
"The area of the semi-circumference: 1.571174; Our approximation of pi: 3.142348; We missed of -0.0240434229862333 percent"

In this case we made an approximation error of -0.02 percent with respect to the built-in value of pi. Not bad huh?

More maths with Python

Today I had some spare time, so I decided to try and write some basic functions in Python which could be useful to have some fun.

Here they are:

Function no. 1 PRIME NUMBERS FINDER

#find prime numbers less than n
def find_primes(n):
    i = 2
    primes = []
    while i <= n:
        if len(find_divisors(i)) == 2:
            primes.append(i)
            i += 1
        else:
            i += 1
    return primes

print(find_primes(1000))

This function finds all the prime numbers which are less than the given number n. Pretty interesting if you are fascinated by prime numbers.

Function no. 2 PRIME FACTORIZATION FUNCTION

#Prime factorization
def factorize(n):
    a = n
    i = 2
    factors = [1]
    while i <= n+1:
        if (a%i == 0):
            factors.append(i)
            a = int(a/i)
            i = 2
        else:
            i += 1
    return factors

print(factorize(90))

This function returns the prime factorization of a given number. As usual, the item returned is a list.

Function no. 3 FIND ALL THE DIVISORS

#Find divisors
def find_divisors(n):
    divisors = [1]
    i = 2
    while i <= n:
        if n % i == 0:
            divisors.append(i)
            i += 1
        else:
            i += 1
    return divisors

print(find_divisors(90))

This function returns the list of the divisors of a given number.

Function no. 4 FIND THE MCD

#find MCD
def find_MCD(a,b):
    divisors_a = find_divisors(a)
    divisors_b = find_divisors(b)
    common_divisors = []
    for i in divisors_a:
        if i in divisors_b:
            common_divisors.append(i)
    return max(common_divisors)

print(find_MCD(20,40))

This function returns the MCD of two given numbers and can be pretty handy if you need to reduce a fraction. Furthermore, it can easily be extended to more arguments.

I wish I could program when I was in high school!! :)

How to permanently store values in python? sqlite3!

When you start programming, the first thing you understand is that there is a ton of things, problems and solutions that you are going to deal with. And, in many cases, you do not need to study from books (though that is still an important part) to get a grasp of a new function or piece of code because you will have a problem to solve and that piece of code is the solution that you need. Perhaps you would not even remember it had it not occurred in such a way.

This is the case for the sqlite3 module and me. When developing small projects in Python, sooner or later you are going to face the problem of how to permanently store a value of a variable. Now, that is fine if you need to store it while the script is running, but as soon as you shut down your pc or mac, the moment you shut down Python, that value goes lost. Unless it is in the script and it will be read again next time you run it, it will disappear. But what if we want to make Python remember data entered from the user for future use?

As far as I know, there are actually a few options:
-An external .txt file
-The module Pickle (though I know nothing else about it so far)
-A database

I shall discuss the third option.
Even though there are different modules for this option, not all of them are available for Python 3 (I refer to mysql).
The module sqlite3 is a ready-to-use, built-in module which provides a set of instruction to make Python interact with a database.

Let's see how to create a database:

import sqlite3

conn = sqlite3.connect("my_database.db")

This piece of code will connect to the database named my_database or create a new one named my_database in case it does not find the database.

Now you need to create an object (I do not know if that is the correct name, in case it is not please correct me) which will help us executing SQL code:

c = conn.cursor()

Every time you want to interact with the database now you need to use c, for instance, to create a table:

c.execute("CREATE TABLE products (id INTEGER, PRIMARY KEY, 
name text, model text, cost int, price real)"), 

The string in the parenthesis is SQL code.
SQL is the language by which you can interact with a database.
If you are already familiar with SQL syntax then you will not have any problem using this module, if you are not then I suggest you to check out SQL tutorials for the basic commands.

To insert data into a table you can just do like this:

c.execute("INSERT INTO products VALUES ('orange','brand1'
,1,2.3)")

And now the trickiest part (for me): query your database and show data
In order to query your database, you could do like this

query1 = c.execute("SELECT * FROM products")
query1

However, this will NOT display any data on your Python command line!! In order to display data, you need to iterate over query1, like this:

for row in query1:
    print(row)

This will display the results of your query and print them out.

After you are done with your database, do not forget to close the connection. Of course, next time you will have to reopen it again.

conn.close()

Hope this was useful.

Maths with Python

Last week I was reviewing analytic geometry and while looking at graphs, equations, theorems and points I was a little bored and decided to try and implement some part of that in Python.

Even though I am still not that practical with Python classes, I am beginning to experiment classes for different purposes. As you know a class is a sort of template, a collection of code and methods (functions inside a class are called methods) which can be called on the instances of the class.

In this example, I created two classes: a class for points and a class for circumferences.
The modules I used are matplotlib and numpy. You can check on my Python page for the link to download these modules or simply google them.

The Point class is really basic, well, not the most basic class you could create but pretty basic.

class Point(object):

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def __repr__(self):
        return "(%s,%s)" %(round(self.x ,2), round(self.y, 2))

The __init__ methods initializes the class. After you have defined an instance of the point Class, as below

p = Point(0,0)

every time you need to know the x coordinate, you can simply type in

p.x

and Python will print to you the value of x for the instance (the point) p.

The __repr__ methods gives you the freedom to decide how to (represent) display your instance when called. For example, in this case, when I type in the command line

p

The interpreter will show (0,0) because that is the representation of a point p I told him to display through __repr__
The screenshots below show the results

Now that we have our Point class, let's skip to fancier stuff :)
for the circumference class pyplot and numpy are needed so we need to import them. Note that usually import statements are above, in the first lines of the script, however, in order not to confuse anyone I put them here.
This is the Circle class:

#New class for circumferences
import matplotlib.pyplot as plt
import numpy as np

class Circle(object):

    def __init__(self, point, r):
        self.point = point
        self.xc = point.x
        self.yc = point.y
        self.r = r
        self.a = -2*self.xc
        self.b = -2*self.yc
        self.c = (self.xc)**2 + (self.yc)**2 - (self.r)**2

    def check_if_circle(self, x, y):
        if (round(x**2 + y**2 + self.a*x + self.b*y + self.c, 1) == 0):
            #print("The point is on the circumference")
            return True
        else:
            #print("The point is NOT on the circumference")
            return False


    def __repr__(self):
        parameters = {"a":self.a, "b":self.b, "c":self.c, "r":self.r}
        #print(parameters)
        if self.a != 0 and self.b != 0:
            return ("The circle has the centre in (%s, %s) and
 has a radius of %d.\nThe equation is:\nx^2 + y^2 " + str(self.a)
 + "x " + str(self.b) + "y " + str(self.c) + " = 0") 
 %(round(self.xc, 2), round(self.yc, 2), self.r)
        elif self.a == 0 and self.b != 0:
            return ("The circle has the centre in (%s, %s) and
 has a radius of %d.\nThe equation is:\nx^2 + y^2 " + str(self.b)
 + "y " + str(self.c) + " = 0")  %(round(self.xc, 2), 
 round(self.yc, 2), self.r)
        elif self.b == 0 and self.a != 0:
            return ("The circle has the centre in (%s, %s) and
 has a radius of %d.\nThe equation is:\nx^2 + y^2 " + str(self.a)
 + "x " + str(self.c) + " = 0")  %(round(self.xc, 2), 
 round(self.yc, 2), self.r)
        else:
            return ("The circle has the centre in (%s, %s) and
 has a radius of %d.\nThe equation is:\nx^2 + y^2 " + str(self.c)
 + " = 0") %(round(self.xc, 2), round(self.yc, 2), self.r)

    def area(self):
        return (np.pi * (self.r)**2)

    def area_sphere(self):
        return (4*(np.pi)*(self.r)**2)

    def volume_sphere(self):
        return(4/3*(np.pi)*(self.r)**3)

    def plot_crf(self):

        def semi_sphere_positive(x):
            y = self.yc + np.sqrt((self.r)**2 - (x - self.xc)**2)
            return y

        def semi_sphere_negative(x):
            y = self.yc - np.sqrt(self.r**2 - (x - self.xc)**2)
            return y

        z = np.linspace(self.xc - self.r,self.xc + self.r,1000)
        graphone = semi_sphere_positive(z)
        graphtwo = semi_sphere_negative(z)
        plt.plot(z,graphone)
        plt.plot(z,graphtwo)
        plt.show()


#define an instance, a circumference of centre p and radius 2
C = Circle(p, 2)

The __repr__ method is slightly fancier in order to display the equation of the circumference in its canonical form according to the values of a,b and c. Remember that:

After this class has been defined, you can call all its methods on C and take a look at the results.
This is particularly interesting if you are actually doing your homework on analytic geometry and you want to develop your programming skills at the same time. Every geometrical figure has its own characteristics and properties. Classes are the perfect tool to collect data which is characterized by a common denominator.

Here you can find a screenshots of the methods in the Circle class

and finally the most interesting method: the graphical representation of the circumference:

C.plot_crf()

and the result:

Note that, since the equation of the circumference is not a function, I had to plot the two semi-circumference. The graphical result looks still good though.
Alternatively, you could have used polar coordinates and the code would have shrunk down by a great amount.

Hello Everyone!

Hi everyone!

My (nick)name is Mic, I am a student and in part of my spare time I enjoy learning programming languanges and practicing my programming knowledge. I have decided to open this blog because no one that I know (or heard of, or that is known to my friends :-)) has the same passion/hobby and therefore I cannot find anyone to share my projects, small programs or views on programming topics with, and last but of course not least, anyone who could point at my programming errors or bad organised code or whatever wrong and clumsy there might be and share their constructive views, opinions and information.

I feel like I have to make some premises:

1. I am not a professional programmer nor a computer science student, I am just a free time programmer who loves technology and sciences, therefore I might as well be doing tons of errors in my code, and most likely not using the appropriate terminology. I apologize for that, I hope this blog and the interaction with the internet community will help me to improve on this point. I might as well not be following the programmer’s best practices, please let me know if you should spot something like this in my code. Hopefully, the quality of my posts will improve over time.

2. The code, pieces of code which I am going to post might have bugs! As the programs I develop grow in the number of functions used and deepness of the code, the more the bugs become difficult to spot and fix. I usually check my code many times before publishing, mainly because I would like it to work as a first step and then perhaps I will publish it. Not every piece of code I write is going to be shared of course.

3. I mostly program in Python, Java and R. I usually use R for data analysis and simulations (as it should be I guess), while Python and Java for more interactive and complex stuff. I also got to learn about C++ whose syntax is similar to the one used in Java.

4. The programs I might post are entirely made by me and by my imagination or simulation of well known games or algorithms. Every reference that might occur is not intended and will be removed if I am asked to do so. I am not responsible for any kind of use people might do with my code. The sole purpose of this blog is educational and the sharing of views, opinions, tips, information among beginner users. Check the disclaimer for more information.

Mic