To begin with, we have two functions. The first is a probability density function (probability mass function in the code implementation, since time must be discretized) which gives the probability of flying at any given time of day assuming that the probabiilty of taking a flight sometime during the day is 1. This function is called f_time_prob in the code, $f(t)$ here in text. In this example we will suppose it looks like this:
In [62]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
Nmin = 1440 # 1440 minutes in a day!
def f_time_prob(n=1440, a1=7,b1=1.5,a2=19,b2=1.5):
"""Returns probability mass function of flying based on time at n equally
spaced times.
Arguments:
- n -- number of wind data points per day available
- a1 -- time of morning at which to return 0.5
- b1 -- steepness of morning curve (larger numbers = steeper)
- a2 -- time of afternoon at which to return 0.5
- b2 -- steepness of afternoon curve (larger numbers = steeper)
Returns:
- A probability mass function for flying in each of n intervals
during the day"""
# time is in hours, and denotes start time of flight.
# wind was recorded starting after the first 30 min, so exclude midnight
# at the start, include it at the end.
t_tild = np.linspace(0+24./n,24,n) #divide 24 hours into n equally spaced times
# Calculate the likelihood of flight at each time of day, giving a number
# between 0 and 1. Combination of two logistic functions.
likelihood = 1.0 / (1. + np.exp(-b1 * (t_tild - a1))) - \
1.0 / (1. + np.exp(-b2 * (t_tild - a2)))
# Scale the likelihood into a proper probability mass function, and return
return likelihood/likelihood.sum()
f_func = f_time_prob()
In [63]:
time_vec = np.linspace(0+24./Nmin,24,Nmin) # time elapsed in hours
plt.figure()
plt.plot(time_vec,f_func*60) # multiplying by 60 here so it looks like the probability density function
# would, matching the x-axis units of hours. In code, f_func is a
# probability mass function - f_func.sum() = 1.
plt.tick_params(axis='both', which='major', labelsize=14)
plt.xlabel('time (hrs)',fontsize=21)
plt.ylabel(r'$f(t)$',fontsize=21)
plt.title('Daylight Probability Density',fontsize=21)
plt.axis([0,24,0,0.09])
plt.show()
In [64]:
def g_wind_prob(windr, aw=1.8,bw=6):
"""Returns probability of flying under given wind conditions.
If the wind has no effect on flight probability, returns 1.
Otherwise, scales down from 1 to 0 as wind speed increases.
Arguments:
- windr -- wind speed
- aw -- wind speed at which flight prob. is scaled by 0.5
- bw -- steepness at which scaling changes (larger numbers = steeper)"""
return 1.0 / (1. + np.exp(bw * (windr - aw)))
# Generate some "wind"
np.random.seed(100)
sigma = 0.1
wind = np.zeros(Nmin)
S = sigma*np.random.randn(Nmin)
wind[0] = np.max((0,S[0]))
for ii in range(1,Nmin):
wind[ii] = np.max((0,wind[ii-1] + S[ii]))
g_func = g_wind_prob(windr=wind)
fig = plt.figure(figsize=(16,4))
plt.subplot(121)
plt.plot(time_vec,wind,'g')
plt.xlabel('time (hrs)',fontsize=18)
plt.ylabel('wind speed (km/hr)',fontsize=18)
plt.title('Wind Speed During the Day',fontsize=20)
plt.axis([0,24,0,4])
plt.subplot(122)
plt.plot(time_vec,g_func)
plt.xlabel('time (hrs)',fontsize=18)
plt.ylabel(r'$g(\mathbf{w}_d)$',fontsize=18)
plt.title('Scaling of Flight Probability by Wind Speed',fontsize=20)
plt.axis([0,24,0,1])
fig.subplots_adjust(wspace=0.3)
plt.show()
To scale $f(t)$ by $g(t,w(t))$, we simply multiply. The result, $f_g(t,w(t))$ is no longer a probability density function (mass function in the descretized code) because it does not integrate/sum to one, but we can simply consider the difference as the probability that no flight is taken during the day.
In [65]:
plt.figure()
plt.plot(time_vec,f_func*60*g_func)
plt.xlabel('time (hrs)',fontsize=18)
plt.ylabel(r'$f(t)\cdot g(\mathbf{w}_d(t))$',fontsize=18)
plt.title(r'$g(\mathbf{w}_d(t))$ applied to $f(t)$',fontsize=20)
plt.show()
In [66]:
# This just for presentation display purposes - combine the first and last g plot
fig, ax1 = plt.subplots(figsize=(12,6))
plt.title('Flight Probability Density',fontsize=24)
ax1.plot(time_vec,f_func*60,color='#FFCC33',alpha=1)
ax1.plot(time_vec,f_func*60*g_func,'k')
ax1.set_xlim(0,24)
ax1.set_ylabel('probability density',fontsize=21)
ax1.set_ylim(0,0.09)
ax1.tick_params(axis='both', which='major', labelsize=14)
ax1.set_xlabel('time (hrs)',fontsize=21)
ax2 = ax1.twinx()
ax2.plot(time_vec,wind,'g',alpha=0.5) # wind speed
ax2.set_xlabel('time (hrs)',fontsize=21)
ax2.set_ylabel('wind speed (km/hr)',color='g',fontsize=21, alpha=0.5)
ax2.set_ylim(0,4)
ax2.tick_params(axis='both', which='major', labelsize=14)
for t2 in ax2.get_yticklabels():
t2.set_color('g')
t2.set_alpha(0.7)
#plt.savefig('flight_prob_density.pdf',format='pdf')
plt.show()
In [71]:
# This cell is to reproduce the above for journal submission. It should be fine both in color and B/W.
alpha_pow = 1
min_vec = np.linspace(1,24*60,24*60)
integral_avg = f_func*g_func/min_vec/np.amax(f_func)*np.cumsum((1-np.cumsum(f_func)**alpha_pow)*(f_func-f_func*g_func))
h_func = f_func*g_func + integral_avg
fig, ax1 = plt.subplots(figsize=(12,6))
plt.title('Wind-borne flight take-off probability density',fontsize=24)
ax1.plot(time_vec,wind,':b',alpha=0.6,label='wind speed') # wind speed
ax1.set_xlabel('time (hrs)',fontsize=18)
ax1.set_ylabel('wind speed (m/sec)',color='b',fontsize=18, alpha=0.6)
ax1.set_ylim(0,4)
ax1.tick_params(axis='both', which='major', labelsize=14)
for t1 in ax1.get_yticklabels():
t1.set_color('b')
t1.set_alpha(0.6)
ax1.legend(loc='upper right',fontsize=14)
ax2 = ax1.twinx()
ax2.plot(time_vec,f_func*60,'--',color='#DDBB22',alpha=1,label='daylight flight probability')
ax2.plot(time_vec,f_func*60*g_func,'c',alpha=0.5,zorder=1,label='scaled probability')
ax2.plot(time_vec,h_func*60,linewidth=2,color='#2fb47c',alpha=0.75,zorder=0,label='redistributed probability')
ax2.legend(loc='upper left',fontsize=14,framealpha=0.8)
ax2.set_xlim(0,24)
ax2.set_ylabel('prob. density',fontsize=18)
ax2.set_ylim(0,0.11)
ax2.tick_params(axis='both', which='major', labelsize=14)
xticks = ax2.get_xticks()
xticks[-1] = 24
ax2.set_xticks(xticks)
ax2.set_xlabel('time (hrs)',fontsize=21)
plt.tight_layout()
#plt.savefig('flight_prob.eps',bbox_inches='tight',format='eps')
plt.show()
The problem with this approach is that in essence, it only gives insects one shot at taking a flight during the day. It's as if each flyer says picks a time in advance according to $f(t)$ and then, once the wind condition at the moment is known, actually takes that flight or not with probability $g(t,w(t))$. If no flight is taken at that moment, none is taken during the day. Modeling flight in this way will most likely underestimate the number of fliers per day significantly.
What we want to do instead is allow those insects who were unable to take a flight early in the day the chance to take a flight later. We will do this in a weighted way - distributing the probability over the rest of the day but not forcing a flight if the day is already late. Another way of looking at it: as the day wears on, there should be less chance that an insect which didn't take a flight earlier will take one now, simply because there is less time left in the day and smaller scale activity may have been taking place earlier. In doing this calculation, we should take into account $f(t)$ and $g(t)$, so that insects are not taking flights at times they normally wouldn't anyway.
The model I have come up for is $$h(t) = f_g(t) + \frac{f(t)g(t)}{f_{max}\cdot t}\int_0^t (1-[\int_0^sf(r)\ dr]^{\alpha})(f(s) - f_g(s))\ ds$$ This calculation is basically a weighted average (where the weight is the cumulative distribution function of $f(t)$), which is then scaled by $g(t)$ and the maximum of $f$. The result is then added back to $f_g(t)$.
In [32]:
# change this value, 0 < cur_s <= 24, to see calculation at different stages.
cur_s = 16
alpha_pow = 1
int_1_minus_f = 1 - np.cumsum(f_func)**alpha_pow
min_vec = np.linspace(1,24*60,24*60)
integral_avg = f_func*g_func/min_vec/np.amax(f_func)*np.cumsum((1-np.cumsum(f_func)**alpha_pow)*(f_func-f_func*g_func))
s_vec = np.linspace(0+24./(cur_s*60),cur_s,cur_s*60)
h_func = f_func*g_func + integral_avg
fig, ax1 = plt.subplots(figsize=(14,8))
ax1.plot(time_vec,f_func*60,'k',label=r'$f(t)$')
ax1.hold(True)
ax1.plot(time_vec,f_func*g_func*60,color='b',label=r'$f_g(t,\mathbf{w}_d)=f(t)\cdot g(\mathbf{w}_d)$')
ax1.fill_between(s_vec,f_func[:s_vec.size]*g_func[:s_vec.size]*60,f_func[:s_vec.size]*60,
facecolor='yellow',alpha=0.5,label='Prob. difference')
ax1.plot(s_vec,h_func[:s_vec.size]*60,color='c',label=r'New probability, $h(t,\mathbf{w}_d)$')
ax1.hold(False)
ax1.set_ylim(0,0.11)
# make sure the y-axis label and tick labels match the line color
ax1.set_xlabel('time (hrs)',fontsize=18)
ax1.set_ylabel(r'$f(t)$, $f(t)\cdot g(\mathbf{w}_d)$, and $h(t,\mathbf{w}_d)$',fontsize=18)
# for t2 in ax2.get_yticklabels():
# t2.set_color('b')
ax1.legend(loc='upper left',fontsize=16)
ax2 = ax1.twinx()
ax2.plot(time_vec,int_1_minus_f,'0.80',label=r'1-$F(t)$')
ax2.set_ylabel(r'$1-F(t)$',fontsize=18)
ax2.set_ylim(0,1)
ax2.hold(True)
ax2.plot(s_vec,int_1_minus_f[:s_vec.size],'m',label=r'1-$F(t)$ weight')
ax2.hold(False)
ax2.legend(loc='upper right',fontsize=16)
plt.title('Redistributing Lost Flight Probability',fontsize=20)
plt.show()
print('Total probability of flying during the day: {0} vs. {1} before.\n'.format(h_func.sum(),(f_func*g_func).sum()))
In [ ]: