Generating Multivariate T Samples in Python

03/12/2013 in Math | Permalink

In order to generate samples from the student’s multivariate t distribution, all you need is some routine to generate multivariate normal samples, and a routine to generate one-dimensional gamma distribution samples. 

In mathematical terms, if you want to sample from a multivariate t distribution t_{N}(\mu,\Sigma) with N degrees of freedom, mean \mu, and scale matrix \Sigma, compute

    \[ \mu + \frac{Z}{\sqrt{\gamma}} \]

where Z \sim \mathcal{N}(0,\Sigma), and \gamma \sim \Gamma(k=N/2,\theta=2/N), where k and \theta are the shape and scale parameters, respectively. 

In python, the code is barely three lines!

import numpy as np
def multivariatet(mu,Sigma,N,M):
    '''
    Output:
    Produce M samples of d-dimensional multivariate t distribution
    Input:
    mu = mean (d dimensional numpy array or scalar)
    Sigma = scale matrix (dxd numpy array)
    N = degrees of freedom
    M = # of samples to produce
    '''
    d = len(Sigma)
    g = np.tile(np.random.gamma(N/2.,2./N,M),(d,1)).T
    Z = np.random.multivariate_normal(np.zeros(d),Sigma,M)
    return mu + Z/np.sqrt(g)

Beginner’s Guide to Git and Github

02/08/2013 in Code | Permalink

Git is a great tool for version control. If you write a little or a lot of code, it’s a good idea to keep track of changes you made. This is extremely useful, especially if you found a mistake and need to undo things. Git also allows you to creat branches which allows you to prototype without worrying about loosing your original, pre-prototyped work. Now, you can do this all without a version control system, but Git allows you to do this all within an organized framework. For example, you could easily just copy your source code to a new folder and try experimenting with a new idea, but if you do this a lot, eventually you will have many different folders, each with different code. In git, you can create branches instead, and only checkout or work on branches one at a time, so your digital workspace isn’t cluttered.

This guide is meant for those who work in the terminal. In particular, this is for anyone working within Linux (e.g. Ubuntu) or Max OS’s where the terminal is available by default. I haven’t used windows in almost a decade, but I’m sure there is a way to access the terminal (I’m not talking about the DOS window). I think you may be able to use Cygwin, but I am not 100 % sure. Anyways, now for the basic basic guide. For a more thorough guide, check out

http://learn.github.com/p/intro.html

which teaches you not only Git, but how to backup your code online with Github and use it thereafter.

Basic guide to Git

Initialize, add and commit

To use git, all you need to do is type in

git init

In the folder you want your version control system. Then, whenever you want to log and keep track of a change, you type in

git add .
git commit

The first command, git add ., stages the changes for git for all files in teh current directory (you can also add individual files or use the * notation). This essentially means that you are ready to commit or log the change. This is not a permanent thing. You can revert any changes at any time (keep reading below). The next command,git commit, actually logs or commits the change. As soon as you type this, the default text editor, e.g. vi, will open up. Write a brief description of the changes made and then save and quit. You can type

git log

To see a few of the last commits. Other useful commands are git status and git diff. You can also make changes to your previous commit by typing

git commit --amend

which is useful if you want to add or remove details about the commit. Also, if you commit without a description, it will not commit.

Read the rest of this entry »

Tags: , ,

Pages: 1 2 3

Animating the planets! Creating animated gifs in python.

01/31/2013 in Math, Physics | Permalink

orbit1

In a nutshell…

Creating animated plots in python is really easy (of course it took me a while to figure it out). We will actually create an animated gif to visualize the motion of the planets. 1 The basic idea is

1) Create a .png file for every frame of your movie using the savefig command from the matplotlib package.

2) Combine all the .png files into an animated gif using ImageMagick.

Step (2) can actually be done within the command line in your terminal, but we will show you a way of doing it all within the python code. You can kip all the details and download it here (compare this code to the un-animated version found here). If you want all the details, keep reading.

The details…

Let’s start with the code introduced in our last post about simulating the motion of the planets. I will insert the code below for convenience.


from numpy import *
from matplotlib.pyplot import *
from scipy.integrate import ode

'''
Numerical solution of the one-body problem using dopri5
'''

# Compute RHS of ODE, f(t,Y)
def f(t,Y):
    x,y,vx,vy = Y # define individual values
    d3 = (x**2+y**2)**1.5 # ||x||^3
    M = .5
    return array([vx,vy,-M*x/d3,-M*y/d3])

# initial parameters
t0 = 0 # initial start time
tfinal = 5 # final start time
dt = .005 # time step to solution (dopri5 algorithm uses adaptive)
y0 = array([.25,0.0,0.0,.45]) # initial position and velocity

# initiate integrator object
test = ode(f).set_integrator('dopri5',atol=1e-6) # prescribe tolerance for adaptive time step
test.set_initial_value(y0,t0) # set initial time and initial value

fig = figure() # initialize figure
ax = fig.add_subplot(111) # name of the plot
while test.successful() and test.t < tfinal:
    # integrate
    test.integrate(test.t+dt)
    # plot position
    ax.plot(test.y[0],test.y[1],'b.',markersize=2)
    ax.plot(0,0,'oy',markersize=12)
    ax.set_xlim([-.15,.35])
    ax.set_ylim([-.12,.15])
    ax.set_title('One-body orbit')
show() # show plot

The code we need to adjust is contained within while loop, highlighted above. We need to save each plot in the while loop to a .png file. This is really easy. The only tricky thing is to create a variable filename (which is not even necessary, but just convenient). I changed the above highlighted code to the following:

Read the rest of this entry »

Notes:

  1. What I like about animated gifs is that they are extremely easy to embed into websites, emails, etc. ↩

Tags: , , , , ,

Pages: 1 2

Simulating the planets with python, Part II: Solving the equations

01/17/2013 in Uncategorized | Permalink

In the previous post, we discussed the theory behind the single-body problem, i.e. of a planet orbiting about a fixed object. Recall that the governing equation is given by

(1)   \begin{equation*} \begin{align} \frac{d^2 \bf{x}}{dt^2} &= G M \frac{\bf{-x}}{\| \bf{x} \|^3}, \end{align} \end{equation*}

where \bf{x} \in \mathbb{R}^2 to start. How do we solve this? One can take the analytic approach, which transforms the equation into a classical central force problem. From here, there are various options, one of which is to solve the resultant ode via a polar coordinate transformation. For our purposes, we will skip the analysis (since we do not want to turn this into an ode exercise) and opt for a numerical approach in python.

We are going to use the ode integration toolbox within scipy to solve (1). Most ode tools work for first-order ode systems, whereas (1) is a second-order ode. What do we do in this case? The standard trick is to turn the second-order ode into a first-order system by introducing the velocity variable. Consider the equivalent formulation of (1) above:

(2)   \begin{equation*} \begin{align} \frac{d \bf{x}}{dt} &= \bf{v}, \\ \frac{d \bf{v}}{dt} &= G M \frac{\bf{-x}}{\| \bf{x} \|^3}, \end{align} \end{equation*}

where we have introduced the velocity \frac{d \bf{x}}{dt} = \bf{v}. Now, (2) is a first-order ode system!

Read the rest of this entry »

Simulating the planets with python, Part I: Theory

01/15/2013 in Uncategorized | Permalink

Have you ever wondered why the orbits are circular or elliptic? Why they are always moving in the same plane? Or perhaps you have seen simulation in movies and wondering why the velocities of the planets speed up as you get closer to the object they are orbiting around. Until now, I have never actually sat down to figure out what the governing equations are and solve them (either analytically or numerically). This post will be about solving a simple one-body problem, i.e. the orbit of a body around a fixed, second body, using a built-in numerical solver in python.

The first step is to determine the governing equations of motion. Let m denote the mass of object 1 and M denote the mass of object 2. For now, let object 2 be fixed. Let \bf{x} and \bf{v} be the position and velocity vectors. For simplicity, we will start with two dimensions so that \bf{x},\bf{v} \in \mathbb{R}^2.  Since object 2 is fixed, the position and velocity components will be both be the zero vector. In summary, we have

Variables: \bf{x},\bf{v} \in \mathbb{R}^2

Constants: m,M

Now, the first piece of the governing equations is Newton’s law of universal gravitation. 1 This law states that the force of one object on the other is inversely proportional (in magnitude) to square of the distance between the two objects. More precisely, the direction of that force is in the direction of the line connecting the two objects (see Figure 1).

For our one body problem, this force on object 1 to the object centered at the origin is then given by

(1)   \begin{equation*} \begin{align} \bf{F} &= G m M \frac{\bf{-x}}{\| \bf{x} \|^3}. \end{align} \end{equation*}

First note that the magnitude of the force is proportional to the inverse square of the distance, i.e., |F| \propto 1/|x|^2 and proportional to the product of their masses . Secondly, note that the direction of the force is in the direction towards object 2 at the origin. Lastly, G is the universal gravitation constant. For simplicity, in our numerical scheme, we will take G,M = 1.0.

Next, Newton’s second law of motion states that this force must satisfy

(2)   \begin{equation*} \begin{align} \bf{F} &= m \times \bf{a}. \end{align} \end{equation*}

where

    \[ \bf{a} = \frac{d^2 \bf{x}}{dt^2}. \]

Equating (1) and (2) gives us our governing equation of motion as a second order, time-dependent, ordinary differential equation, in two dimensions:

(3)   \begin{equation*} \begin{align} \frac{d^2 \bf{x}}{dt^2} &= G M \frac{\bf{-x}}{\| \bf{x} \|^3}. \end{align} \end{equation*}

And that’s it! This simple equation will give us the motion of a planet about a fixed mass. Of course, we in order to close this ode, we need two initial conditions. In our next post, we will solve this via a numerical method and simulate our first planetary orbit.

Notes:

  1. See https://en.wikipedia.org/wiki/Newton%27s_law_of_universal_gravitation  ↩

« Older entries