```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import itertools

plt.rcParams["figure.figsize"] = (15, 15)

%matplotlib inline
```

# Modeling the random walk  

We assume the following:

1) X and Y are independent

2) Each jump of the fire over given timestep t

Our approach is to treat the spread of fire as a random walk on $x$ and $y$.
For each step of the random walk, we'll treat the velocity as a random variable from a distribution $N(\mu, \sigma^2)$

Our input data is a set of trajectories based off of the given tweets. 
Each trajectory gives us a set of jumps associated with the given fire of interest in our region. 



```python
numtweets = 30
r = 50
variance = 5
binsize = 4

#x, y, and timestamp of each tweet
testx = np.linspace(0, r, numtweets) + np.random.normal(0,variance,numtweets)
testy = np.linspace(0, r, numtweets) + np.random.normal(0,variance,numtweets)
testt = np.arange(testx.shape[0])


bins = np.arange(min(testt), max(testt), binsize) #make some bins based on binsize
inds = np.digitize(testt, bins) #for each time in testt, assign it to a bin; 
#the bin of element i in testt will be the bin for element i in testx and testy also
 
numbins = len(bins)

x_trajectories = [] 
y_trajectories = [] 

for i in range(numbins):
    tempx = []
    tempy = []
    for j in range(len(testt)):
        if inds[j] == i + 1: #since bin 0 is indexed as 1
            tempx += [testx[j]] #if the jth element of testt is in bin i, 
            tempy += [testy[j]] #add testx[j] and testy[j] to their appropriate bins
    x_trajectories += [tempx]
    y_trajectories += [tempy]

    
combinatorial_x_trajectories = list(itertools.product(*x_trajectories))
combinatorial_y_trajectories = list(itertools.product(*y_trajectories))

x_distances = [np.diff(i) for i in combinatorial_x_trajectories]
y_distances = [np.diff(i) for i in combinatorial_y_trajectories]

x_distances = np.reshape(x_distances, -1)
y_distances = np.reshape(y_distances, -1)

vx = x_distances/binsize
vy = y_distances/binsize
```


```python
def MLE_mu(data):
    n = len(data)
    return 1/n * np.sum(data)

def MLE_sigma2(data):
    n = len(data)
    mu_hat = MLE_mu(data)
    return 1/n * np.sum((data - mu_hat)**2)

```

Since the behavior of fire spread is uncertain, we assume each trajectory represents an equally likely path of the the fire.
Based off of this, each link $A\to B$ gives us a "representative" sample of the underlying distribution  of velocities, both of x and y.

Therefore, the approach will to be to calculate the $v_x = d_x/t$ for each link of Tweets per trajectory, then to use these to calculate MLEs for the normal distribution governing the velocities of each hop. 

Once we have the normal distribution for these hops, we can use this to predict the probability that the fire has reached some point $A$ by some time $T$.
There are two ways we can do this:

1) Set some timestep $t$ within the range of timesteps that we have seen in the data, and subdivide the desired segment into even chunks such that $\sum t_i = T$ and then add the normals such that the add up to the correct value. 
In this case, say we have $n$ chunks, 
Then, the probability that a fire has spread to at least to point $A$ is 
$$
\begin{align}
P\left(\sum_{i=1}^n x_i \geq A \right) &= P\left(\sum_{i=1}^n tv_i \geq A \right) \\
&= P\left(\sum_{i=1}^n v_i \geq \frac{A}{t} \right)\\
&= P\left(N(n\mu, n\sigma^2)\geq \frac{A}{t} \right)\\
&= P\left(N(n\mu, n\sigma^2)\geq \frac{A}{T/n} \right)\\
&= P\left(N(n\mu, n\sigma^2)\geq \frac{nA}{T} \right)\\
&= P\left(N(\mu, \sigma^2)\geq \frac{A}{T} \right)
\end{align}
$$

2) Find the average velocity required to traverse the whole path in one go and find that probability.
$$
\begin{align}
P(X \geq A) &= P\left(vT \geq A\right)\\
&= P\left(v \geq \frac{A}{T}\right)\\
&= P\left(N(\mu, \sigma^2) \geq \frac{A}{T}\right)\\
\end{align}
$$

Let's apply these ideas below.


First, calculate velocity components based off of the given data:


```python
vx, vy
```




    (array([-0.06722391,  4.897783  , -0.53739367, ...,  4.91297291,
             2.56377293, -1.20113813]),
     array([ 2.65637064,  0.955601  ,  1.38734809, ..., -0.7181989 ,
             3.03220081, -1.14865946]))




```python
timesteps = np.diff(testt)
vx = np.diff(testx)/timesteps
vy = np.diff(testy)/timesteps
```

Get MLE estimates for $v_x$, $v_y$


```python
muhat_y, sigmahat_y = MLE_mu(vy), MLE_sigma2(vy)
muhat_x, sigmahat_x = MLE_mu(vx), MLE_sigma2(vx)
print("means, x and y: " , muhat_x, muhat_y)
print("variances, x and y: ", sigmahat_x, sigmahat_y)
```

    means, x and y:  1.749640387263358 1.3795732118077177
    variances, x and y:  64.0309180225004 45.88646069867259



```python
vy_dist = norm(loc = muhat_y,  scale = np.sqrt(sigmahat_y))
vx_dist = norm(loc = muhat_x, scale = np.sqrt(sigmahat_x))
```


```python
predictx = [vx_dist.mean() * t for t in testt]
```


```python
plt.scatter(testt, testx, label = "x test")
plt.scatter(testt, predictx, label = "x predict")
# plt.scatter(testt, testy, label = "y test")
plt.legend();
```


```python
trajectories = [
    [t * vx_dist.rvs() for t in testt]
    for i in range(5)
]

for t in trajectories:
    plt.plot(testt, t)
```


```python
simulated_y_trajectories = []
for i in range(40):
    simulated_position = [0]
    last = 0

    for t in testt[1:]:
        last += vy_dist.rvs()
        simulated_position += [last]
    
    plt.plot(testt, simulated_position)
    simulated_y_trajectories += [simulated_position]
```


```python
simulated_x_trajectories = []
for i in range(40):
    simulated_position = [0]
    last = 0

    for t in testt[1:]:
        last += vx_dist.rvs()
        simulated_position += [last]
    
    plt.plot(testt, simulated_position)
    simulated_x_trajectories += [simulated_position]
```


```python
import seaborn as sns
simulated_x_trajectories_plot = np.reshape(simulated_x_trajectories, -1)
simulated_y_trajectories_plot = np.reshape(simulated_y_trajectories, -1)
```


```python
plt.figure(figsize = (10, 10))
sns.kdeplot(simulated_x_trajectories_plot, simulated_y_trajectories_plot, shade = True)
for x_traj, y_traj in zip(simulated_x_trajectories, simulated_y_trajectories):
    plt.plot(x_traj, y_traj)
plt.scatter(testx, testy, label = "test data", color = 'orange')
plt.legend();
```

# What happens with multiple fires?

Inclusion exclusion:
$$
\begin{align}
P(fire @  A @  t) &= P(fire 1 @  A @  t) + P(fire 2 @  A @  t) - P(fire 1 \cap fire 2 @  A @  t)\\
&= P(fire 1 @  A @  t) + P(fire 2 @  A @  t) - P(fire 1 @  A @  t)P(fire 1 @  A @  t)
\end{align}
$$
for time $t$ and point of interest $A$


```python

```

c0080f27e5051f9abe33569c373bf4609a1295f4
ipynb
Jupyter Notebook
randomwalk.ipynb
SyedAbutalib/Team-X
e103ac733732782bb3e5cdfc0b578b53ad179ad2
MIT
2019-11-02T20:11:38.000Z
2019-11-03T00:22:29.000Z
randomwalk.ipynb
SyedAbutalib/FireWatch
e103ac733732782bb3e5cdfc0b578b53ad179ad2
MIT
randomwalk.ipynb
SyedAbutalib/FireWatch
e103ac733732782bb3e5cdfc0b578b53ad179ad2
MIT
2019-11-04T13:35:35.000Z
2019-11-04T13:35:35.000Z
Qwen/Qwen-72B
1. YES
2. YES