\documentclass[10pt,vi,fullpage,graphicx,letter]{article}
\usepackage{fullpage}
\usepackage{graphicx}
\title{Nature of Mathematical Modeling\\
PS2}
\author{Piotr Mitros}
\begin{document}
\maketitle
Worked with Andy.

\noindent
7.1.a) I don't understand this question. By definition, the $n$th cumulant is: 
$$C_n=\frac{d^n \Psi(0)}{ds^n}$$
Where $\Psi$ is the second characteristic function. They can also be
expressed as:
$$C_1=\mu_1=E(X)$$
$$C_2=\sigma^2$$
$$C_3=\mu_3=E(x^3)$$

\noindent
7.1.b) 

$$\phi(s)=e^{\mu s +\sigma^2 s^2/2}$$
$$\psi(s)=\mu s +\frac{\sigma^2 s^2}{2}$$

$$C_0= \mu s +\frac{\sigma^2 s^2}{2} = 0$$
$$C_1=\mu +\sigma^2 s = \mu$$
$$C_2=\sigma^2$$
$$C_3=0$$

\noindent
7.2.a) It is: 

$$
\left|
\frac{\partial y_1}{\partial x_2}(x_2) dx_2
\frac{\partial y_2}{\partial x_1}(x_1) dx_1
- 
\frac{\partial y_2}{\partial x_2}(x_1,x_2) dx_2
\frac{\partial y_1}{\partial x_1}(x_1,x_2) dx_1
\right|
$$

So it is scaled by: 

$$
\left|
\frac{\partial y_1}{\partial x_2}(x_1,x_2)
\frac{\partial y_2}{\partial x_1}(x_1,x_2)
- 
\frac{\partial y_2}{\partial x_2}(x_1,x_2)
\frac{\partial y_1}{\partial x_1}(x_1,x_2)
\right|
$$

\noindent
7.2.b) 

$$
\left|
\frac{\partial y_1}{\partial x_2}(x_1,x_2)
\frac{\partial y_2}{\partial x_1}(x_1,x_2)
- 
\frac{\partial y_2}{\partial x_2}(x_1,x_2)
\frac{\partial y_1}{\partial x_1}(x_1,x_2)
\right|=
$$
$$
\left|
\sqrt{-2 \log x_1} \cos (x_2)
\frac{-1}{x_1\sqrt {-2 \log x_1}} \cos (x_2)
+
\sqrt{-2 \log x_1} \sin (x_2)
\frac{-1}{x_1\sqrt {-2 \log x_1}} \sin (x_2)
\right|=
$$
$$
\left|
\frac{-\cos^2 (x_2)}{x_1} 
+
\frac{-\sin^2 (x_2)}{x_1}
\right|=
$$
$$
\frac{1}{x_1} 
$$

Now, we have:
$$y_1^2+y_2^2=-2 \log x_1$$
$$x_1=e^{-\frac{y_1^2+y_2^2}{2}}$$

Since we're scaling inversely to the size of a differential element,
the new distribution is: 

$$p(y_1,y_2)=\left\{
\begin{array}{ccc}
C n e^{-\frac{y_1^2+y_2^2}{2}}  & : & y_1^2+y_2^2 > -2 \log (x_{2,\max}) \\
0 & : & y_1^2+y_2^2 \leq -2 \log (x_{2,\max})
\end{array}
\right.
$$

Where $C$ is a normalization constant, $n$ is the number of times that
the initial distribution is mapped to the new distribution (since
$\sin(x_2)$ and $\cos(x_2)$ are periodic). If $x_1 \in (0,1)$ and $x_2
\in n [0,2\pi) + z$, then this becomes a standard 2d multivariate
normal distribution. 

\noindent
7.2.c) Neil said to just do this in one dimension. Here, 
$$C_1 \approx -0.0013349302731746232 \approx 0$$
$$C_2 \approx 0.9996056556376719 \approx 1$$
$$C_3 \approx -0.0082527053120206185 \approx 0$$

\begin{verbatim}import sys, os, os.path
from pylab import *

seed=262093818;


def randbit():
    global seed
    seed=seed<<1
    if(seed&(2**2)):
        seed=seed^1;
    if(seed&(2**29)):
        seed=seed^1;
    seed=seed&(2**30-1)
    return seed&1;

def randnum():
    sum=0
    for i in range(0,30):
        sum=sum+randbit()
        sum=sum<<1
    sum=sum+0.0
    sum=sum/(2**31)
    return sum;

x1=range(0,100000);
x2=range(0,100000);
y1=range(0,100000);
y2=range(0,100000);

mu1=0;
mu2=0;
mu3=0;


for i in range(1,100000):
    x1[i]=randnum()+randnum()/(2**32)
    x2[i]=randnum()*2*pi
    y1[i]=sqrt(-2*log(x1[i]))*sin(x2[i])
    y2[i]=sqrt(-2*log(x1[i]))*cos(x2[i])
    mu1=mu1+y1[i];
    mu2=mu2+(y1[i]+0.0013349302731746232)**2;
    mu3=mu3+y1[i]**3;

plot(y1, y2, "r+");
show()\end{verbatim}

\noindent
7.3.a) $x_i=x_{i-1}+x_{i-4}$ gives 100011110101100. Generated by: 
\begin{verbatim}  int c=1;
  while(1) {
    c<<=1;
    if(c&2) c^=1;
    if(c&16) c^=1;
    printf("%i", c&1);
  }\end{verbatim}

\noindent
7.3.b) 
$$2^b-1=10^9 \cdot 60 \cdot 60 \cdot 24 \cdot 365 \cdot 10^{10}$$
$$b=9\cdot\log_2(10)+2\cdot\log_2(60)+\log_2(24)+10\cdot\log_2(10)+\log_2(365)$$
So $b>88.027$ and so the shortest one is 89 bits. 

\noindent
7.4.a) 
$$\frac{\delta p}{\delta t}=D\frac{\delta^2 p}{\delta x^2}$$

$$\frac{\delta P}{\delta t}=DPs^2$$
$$\frac{1}{P}\delta P=Ds^2 \delta t$$
Integrating:
$$\log P=Dt s^2+C$$
$$P=e^{Dt s^2+C}$$
$$P=Ce^{Dt s^2}$$
By initial conditions,
$$P=e^{Dt s^2}$$
By the inverse Fourier Transform,
$$p(t,x)=\frac{1}{\sqrt{4 \pi Dt}} e^{-\frac{x^2}{4Dt}}$$

\noindent
7.4.b) It is Gaussian, so $2Dt$

\noindent
7.4.c) The variance goes as $2Dt$ in one form, and as $\frac{kT}{3 \pi
    \mu a} t$ in the other. Therefore,
$$D=\frac{kT}{6 \pi \mu a}$$
Where $\mu$ is the viscosity.

\noindent
7.4.d) 
\begin{verbatim}import sys, os, os.path
from pylab import *

seed=28995;

def randbit():
    global seed
    seed=seed<<1
    if(seed&2):
        seed=seed^1;
    if(seed&32768):
        seed=seed^1;
    return seed&1;

x=range(0,999)
y=range(0,999)
y[0]=0

for j in range(1,10):
    for i in range(1,999):
        y[i]=y[i-1]+randbit()*2-1
    plot(x,y)

for i in range(1,999):
    y[i]=sqrt(i)*3

plot(x,y)

for i in range(1,999):
    y[i]=-sqrt(i)*3

plot(x,y)

show()\end{verbatim}

\includegraphics[scale=0.6]{image.ps}

\noindent
7.4.e) 

It will be: 
$$\frac{1}{1000}
\sum_{t=1}^{1000}
\sum_{i=\left\lceil\frac{t}{2}-
 3 \sqrt{0.25 t}\right\rceil
}^{\left\lfloor \frac{t}{2}+
 3 \sqrt{0.25 t}\right\rfloor
}
\left(
\begin{array}{c}
t \\
i
\end{array}
\right)
0.5^t=
1-
\frac{2}{1000}
\sum_{t=1}^{1000}
\sum_{i=\left\lceil \frac{t}{2} + 3 \sqrt{0.25 t}\right\rceil}^{t}
\left(
\begin{array}{c}
t \\
i
\end{array}
\right)
0.5^t=99.7409
$$ 

\begin{verbatim}warning off MATLAB:nchoosek:LargeCoefficent
sum=0
for t=1:1000
  for i=ceil(t/2-3*sqrt(0.25*t)):floor(t/2+3*sqrt(0.25*t))
    if (i>=0) && (i<=t)
      sum=sum+nchoosek(t,i)*0.5^t;
    end
  end
end
sum\end{verbatim}

The continuous approximation would give 99.73. 

\end{document}