Generate a gaussian wake



In [1]:

    
%pylab inline
import pymc as mc









    



Populating the interactive namespace from numpy and matplotlib

Generate the data



In [2]:

    
import scipy.stats as stats
import numpy as np
figsize(12.5,4)

# Initial sampling space of wind direction [deg]
#x0 = np.linspace(-40,40,10000)
x0 = stats.uniform(loc=-40.0, scale=80.0).rvs(10000)

def wake(teta, teta_e, du_e):
    """ Function defining the wake deficit """
    ## The wake width
    sigma = 30.0
    ## The wake center
    loc = 0.0
    ## The wake deficit
    du = 0.47
    return (du+du_e())*np.exp(-(teta+teta_e()-loc)**2.0/(2*sigma))

wd_uncertainty = 5
ws_uncertainty = 0.00

y0 = array(map(lambda y: wake(y, stats.norm(loc=0.0, scale=wd_uncertainty).rvs, stats.norm(loc=0.0, scale=ws_uncertainty).rvs), x0))
y0_nonoise = array(map(lambda y: wake(y, lambda:0, lambda:0), x0))

# Do a sliding windows averaging
x1 = np.linspace(-40,40,41)
du = x1[1]-x1[0]
y1 = map(lambda x: y0[find(abs(x0-x)<du)].mean(), x1)
# Do the sliding window averaging over the original gaussian
y1_nonoise = map(lambda x: y0_nonoise[find(abs(x0-x)<du)].mean(), x1)

## Plotting
plot(x0, y0, 'k.', label='original data ($\Delta\\theta=%3.1f^\circ, \Delta u/u_\infty =%3.2f$)'%(wd_uncertainty, ws_uncertainty))
plot(x0, y0_nonoise, 'r.', linewidth=2, label='original, without noise')
plot(x1, y1, 'co', markersize=10, label='Sliding window $\\theta=%3.1f^\circ$'%du)
plot(x1, y1_nonoise, 'bo', markersize=10, label='Sliding window (without noise) $\\theta=%3.1f^\circ$'%du)

legend(loc=0)
ylabel('Power Deficit $1-\\frac{P_2}{P_1}$ [-]')
xlabel('Relative wind direction [$^\circ$]')









    Out[2]:





<matplotlib.text.Text at 0x10c59d290>

Help!



In [11]:

    
plot?

Autocomplete



In [ ]:

    
y0_non

Introspection



In [ ]:

    
x0.

Make a horizontal distribution of the power ratio

Plot the wind direction probability density function given a power deficit ratio



In [4]:

    
power_ratios = linspace(0, 0.6, 7)
dpr = 0.1

### Function to find the points within a power ratio +/- power ratio delta
pdf_x0 = lambda prat, pres: x0[find(abs(y0 - prat) < pres)]
for power_ratio in power_ratios:
    try:    
        p = hist(pdf_x0(power_ratio, dpr), 50, normed=True, label='dpr='+str(power_ratio), histtype='step')
    except:
        pass
legend()
ylabel('Pdf')
xlabel('Relative wind direction [$\circ$]')









    



/Users/pire/anaconda/envs/py27/lib/python2.7/site-packages/matplotlib/axes.py:8037: UserWarning: 2D hist input should be nsamples x nvariables;
 this looks transposed (shape is 0 x 1)
  'this looks transposed (shape is %d x %d)' % x.shape[::-1])






    Out[4]:





<matplotlib.text.Text at 0x10c759cd0>



In [5]:

    
hist?

Plot the pdf of of power deficit ratio given a wind direction bin



In [6]:

    
wind_dir = 0.0
bin_sizes = [1, 2.5, 5, 7.5, 10, 15]
pdf_y0 = lambda x, pres: y0[find(abs(x0 - x) < pres)]
for bin_size in bin_sizes:
    try: 
        pdf, bins, patches = hist(pdf_y0(wind_dir, bin_size), 50, normed=True, histtype='step', label='$\Delta^\circ$=%2.1f'%(bin_size))
    except:
        pass
ylabel('Pdf')
xlabel('Power deficit ratio $P_2$/$P_1$ [-]')
legend()









    Out[6]:





<matplotlib.legend.Legend at 0x10ce515d0>

Calculate the bounds of the data using the cumulative probability distribution:

So we want to identify where is located 85% of the data



In [7]:

    
prob, nb =  histogram(pdf_y0(wind_dir, 2.5), 50)

# Cumulative function: Add up the probability
cum = lambda x: array([sum(x[:i+1])/sum(1.0*x) for i in range(len(x))])
# Find where the cumulative function is equal to the bound, given the probability
bounds = lambda x, n, pres: n[argmin(abs(cum(x) - pres))]
# Find the bounds given the wind direction and 
integratedbounds = lambda wd, wdbin, n, pres: bounds(histogram(pdf_y0(wd, wdbin), n)[0], histogram(pdf_y0(wd, wdbin), n)[1], pres) 
    
bound = 0.85
plot(bounds(prob, nb, bound)*ones([2]), [0, 0.1])
plot(bounds(prob, nb, 1-bound)*ones([2]), [0, 0.1])
plot(nb[:-1], prob/sum(1.0*prob), fillstyle='bottom')
indexes = find((nb >= bounds(prob, nb, 1-bound))  & (nb <= bounds(prob, nb, bound)))
fill_between(nb[indexes], 0, prob[indexes]/sum(1.0*prob), color="red", alpha = 0.4 )









    Out[7]:





<matplotlib.collections.PolyCollection at 0x10cf879d0>

Generating the envelops from the artificial data



In [8]:

    
bound = 0.95
x3 = linspace(-30, 30, 100)
plot(x0, y0, 'k.', label='original data', alpha = 0.3)
fill_between(x3, map(lambda wd: integratedbounds(wd, 2.5, 100, 0.75), x3), 
                 map(lambda wd: integratedbounds(wd, 2.5, 100, 1-0.75), x3), 
             color="blue", alpha = 0.4, label='Confidence %f'%bound)
fill_between(x3, map(lambda wd: integratedbounds(wd, 2.5, 100, 0.95), x3), 
                 map(lambda wd: integratedbounds(wd, 2.5, 100, 1-0.95), x3), 
             color="green", alpha = 0.4, label='Confidence %f'%bound)

plot(x0, y0_nonoise, 'r.', linewidth=2, label='original gaussian')
legend()
ylabel('Power Deficit $1-\\frac{P_2}{P_1}$ [-]')
xlabel('Relative wind direction [$^\circ$]')









    Out[8]:





<matplotlib.text.Text at 0x10cf874d0>

Make a boxplot of the datasets



In [9]:

    
# Do a sliding windows averaging
x1 = np.linspace(-40,40,41)
du = 2
y1 = map(lambda x: y0[find(abs(x0-x)<du)].mean(), x1)
plot(x0, map(lambda y: wake(y, lambda:0, lambda:0), x0), 'k.', linewidth=3, label='original, without noise')
plot(x1, y1, 'co-', linewidth=3, label='Sliding window $\\theta$=%3.1f$^\circ$'%(2*du))
b = boxplot(array(map(lambda x: y0[find(abs(x0-x)<du)], x1)).T, positions=x1)
# Plot the confidence interval
bound = 0.75
x3 = linspace(-40, 40, 100)
#fill_between(x3, map(lambda wd: integratedbounds(wd, 2*du, 100, bound), x3), map(lambda wd: integratedbounds(wd, 2*du, 100, 1-bound), x3), color="gray", alpha = 0.4, label='Confidence %f'%bound)
legend(loc=0)
xlim([-40,40])
ylabel('Power Deficit $1-\\frac{P_2}{P_1}$ [-]')
xlabel('Relative wind direction [$^\circ$]')









    Out[9]:





<matplotlib.text.Text at 0x10ce49350>



In [ ]: