This notebook has my attempts at applying Bayesian parameter estimation onto the problem of source-level estimation. 


```python
import arviz as az
import matplotlib.pyplot as plt
import numpy as np 
np.random.seed(100)
import pymc3 as pm
import pandas as pda
from scipy.special import jv as bessel_firstkind
from sympy import besselj, symbols, hankel2, legendre, sin, cos, summation, I, oo, diff, pi
from sympy import Abs, lambdify
```


```python
theta_values = np.linspace(0,np.pi/2,100)
v_sound = 330 # m/s
freq = 50*10**3 # Hz
k_value = 2*np.pi/(v_sound/freq)
a_value = 5*10**-3

# choose another higher frequency
freq2 = 80*10**3
k_value2 =2*np.pi/(v_sound/freq2)

```


```python
# from the online repo: https://github.com/thejasvibr/bat_beamshapes/

def absdivide_or_zero(a,b):
    try:
        out = abs(a/b)
    except:
        out = 0
    return out

def piston_in_infinite_baffle(theta, k=k_value, a=a_value):
    '''
    Generate predictions for the classic piston vibrating in an infinite baffle.

    Parameters
    ----------
    theta: 0<=float.<=2*pi
           The angle in azimuth/elevation at which the relative sound pressure is calculated.
           The relative soundpressure is calculated in comparison to the sound pressure level
           on-axis (0degrees/radians)
    k : float>0,optional
        The wavenumber of the sound (2*pi/wavelength)
    a : float>0, optional
        The diameter/radius(?) of the piston in m. 
    Returns
    --------
    rp_theta : float >=0.
               The relative sound pressure level at the given theta
    '''
    ka_sintheta = k*a*np.sin(theta)
    numerator = 2*bessel_firstkind(1, ka_sintheta)#2*bessel_firstkind(1, ka_sintheta)   
    try:
        rp_theta = np.apply_along_axis(absdivide_or_zero,0,numerator,ka_sintheta)        
    except:
        rp_theta = absdivide_or_zero(numerator,ka_sintheta)
    return rp_theta


#PIBF = lambdify(piston_in_infinite_baffle) 

db = lambda X: 20*np.log10(X)
```


```python

```


```python

rel_levels = np.zeros((theta_values.size,3))
rel_levels[:,0] = theta_values
rel_levels[:,1] = piston_in_infinite_baffle(theta_values)
rel_levels[:,2] = db(rel_levels[:,1])

```

    <ipython-input-40-2b7b0e36b7d7>:5: RuntimeWarning: invalid value encountered in true_divide
      out = abs(a/b)



```python
%matplotlib inline
```


```python
fig = plt.figure()
ax = plt.subplot(111, projection='polar')
ax.set_theta_zero_location("N")  # theta=0 at the top
plt.plot(rel_levels[:,0],rel_levels[:,2])
```

This is a coarse idea of the how beam-shape affects the calculated source-level of bat calls. Now, assuming this the beam-shape in general, let us sample a set of random-points, and generate a set of source-level measurements with beam-shape effects included. 


```python
on_axis_SL = 100 # dB SPL re 20muPa at 1m -   the mean source-level
SL_sd = 1 # a +/- 1dB call-to-call standard deviation

random_thetas = np.random.choice(np.arange(rel_levels.shape[0]),25)
# generaate 
measured_SL = np.random.normal(on_axis_SL,SL_sd,random_thetas.size) + rel_levels[random_thetas,2]
```


```python
plt.figure()
plt.hist(measured_SL);
```


```python
az.summary(measured_SL,kind='stats')
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_3%</th>
      <th>hdi_97%</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>x</th>
      <td>87.99</td>
      <td>8.701</td>
      <td>75.0</td>
      <td>99.796</td>
    </tr>
  </tbody>
</table>
</div>




```python
plt.figure()
az.plot_kde(measured_SL)
```

#### Beam-shape effects can have pretty strong results. 
As we can see above, the beam-shape effects could mean that the measured source-level will be generally under-estimated. The mean is ~82 dB SPL, while the range is between 60-100 dB SPL! 

#### How do we estimate the underlying mean source-level given all the prior information at hand?
$Measured\:level = On-axis\:source\:level - beam\:shape\:effect$ 

Here, the two parts are described by:

* $On-axis\:source\:level = N(\mu,\sigma)$
* $beam\:shape \:effect = f(\theta)$


```python
def fn_measured_level(on_axis,theta):
    return on_axis-db(piston_in_infinite_baffle(theta))
```


```python
m00 = pm.Model()
with m00:
    theta = pm.Uniform('theta',lower=0,upper=np.pi/2,shape=measured_SL.size)
    on_axis = pm.Normal('on_axis',mu=100,sigma=1,shape=measured_SL.size)
    sim = pm.Simulator('sim',fn_measured_level,params=(on_axis,theta),
                                   observed=measured_SL, epsilon=10)
    # -- delete the parallel=True as it always caused issues initially
    trace_lv = pm.sample_smc(draws=500,kernel="ABC")
    idata_lv = az.from_pymc3(trace_lv)
    
```

    Initializing SMC sampler...
    Sampling 2 chains in 2 jobs
    Stage:   0 Beta: 0.052
    Stage:   1 Beta: 0.112
    Stage:   2 Beta: 0.176
    Stage:   3 Beta: 0.246
    Stage:   4 Beta: 0.321
    Stage:   5 Beta: 0.405
    Stage:   6 Beta: 0.496
    Stage:   7 Beta: 0.621
    Stage:   8 Beta: 0.783
    Stage:   9 Beta: 0.987
    Stage:  10 Beta: 1.000
    Stage:   0 Beta: 0.056
    Stage:   1 Beta: 0.120
    Stage:   2 Beta: 0.187
    Stage:   3 Beta: 0.255
    Stage:   4 Beta: 0.327
    Stage:   5 Beta: 0.420
    Stage:   6 Beta: 0.519
    Stage:   7 Beta: 0.640
    Stage:   8 Beta: 0.787
    Stage:   9 Beta: 0.984
    Stage:  10 Beta: 1.000



```python
az.plot_trace(idata_lv, kind="rank_vlines");
```


```python
out_data = idata_lv.to_dataframe()
```


```python
r = az.summary(trace_lv,kind='stats',hdi_prob=0.89)
```

    /home/autumn/anaconda3/envs/statsre/lib/python3.8/site-packages/arviz/data/io_pymc3.py:88: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
      warnings.warn(



```python
theta_summary = r.loc['theta[0]':,:]
theta_summary['hdi_width'] = theta_summary['hdi_94.5%']-theta_summary['hdi_5.5%']
theta_summary['hdiwidth_deg'] = np.degrees(theta_summary['hdi_width'])
theta_summary['mean_deg'] = np.degrees(theta_summary['mean'])
theta_summary['actual_deg'] = np.degrees(theta_values[random_thetas])

theta_summary
```

    <ipython-input-58-6461f91ae797>:2: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      theta_summary['hdi_width'] = theta_summary['hdi_94.5%']-theta_summary['hdi_5.5%']
    <ipython-input-58-6461f91ae797>:3: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      theta_summary['hdiwidth_deg'] = np.degrees(theta_summary['hdi_width'])
    <ipython-input-58-6461f91ae797>:4: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      theta_summary['mean_deg'] = np.degrees(theta_summary['mean'])
    <ipython-input-58-6461f91ae797>:5: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      theta_summary['actual_deg'] = np.degrees(theta_values[random_thetas])





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_5.5%</th>
      <th>hdi_94.5%</th>
      <th>hdi_width</th>
      <th>hdiwidth_deg</th>
      <th>mean_deg</th>
      <th>actual_deg</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>theta[0]</th>
      <td>0.235</td>
      <td>0.185</td>
      <td>0.010</td>
      <td>0.473</td>
      <td>0.463</td>
      <td>26.527946</td>
      <td>13.464508</td>
      <td>40.000000</td>
    </tr>
    <tr>
      <th>theta[1]</th>
      <td>0.253</td>
      <td>0.194</td>
      <td>0.007</td>
      <td>0.501</td>
      <td>0.494</td>
      <td>28.304115</td>
      <td>14.495832</td>
      <td>36.363636</td>
    </tr>
    <tr>
      <th>theta[2]</th>
      <td>0.216</td>
      <td>0.141</td>
      <td>0.029</td>
      <td>0.425</td>
      <td>0.396</td>
      <td>22.689129</td>
      <td>12.375888</td>
      <td>60.909091</td>
    </tr>
    <tr>
      <th>theta[3]</th>
      <td>0.213</td>
      <td>0.157</td>
      <td>0.009</td>
      <td>0.429</td>
      <td>0.420</td>
      <td>24.064227</td>
      <td>12.204001</td>
      <td>65.454545</td>
    </tr>
    <tr>
      <th>theta[4]</th>
      <td>0.273</td>
      <td>0.191</td>
      <td>0.014</td>
      <td>0.545</td>
      <td>0.531</td>
      <td>30.424059</td>
      <td>15.641748</td>
      <td>32.727273</td>
    </tr>
    <tr>
      <th>theta[5]</th>
      <td>0.444</td>
      <td>0.341</td>
      <td>0.008</td>
      <td>0.903</td>
      <td>0.895</td>
      <td>51.279723</td>
      <td>25.439326</td>
      <td>3.636364</td>
    </tr>
    <tr>
      <th>theta[6]</th>
      <td>0.297</td>
      <td>0.221</td>
      <td>0.015</td>
      <td>0.561</td>
      <td>0.546</td>
      <td>31.283496</td>
      <td>17.016847</td>
      <td>37.272727</td>
    </tr>
    <tr>
      <th>theta[7]</th>
      <td>0.366</td>
      <td>0.269</td>
      <td>0.014</td>
      <td>0.705</td>
      <td>0.691</td>
      <td>39.591384</td>
      <td>20.970255</td>
      <td>11.818182</td>
    </tr>
    <tr>
      <th>theta[8]</th>
      <td>0.303</td>
      <td>0.217</td>
      <td>0.019</td>
      <td>0.569</td>
      <td>0.550</td>
      <td>31.512679</td>
      <td>17.360621</td>
      <td>22.727273</td>
    </tr>
    <tr>
      <th>theta[9]</th>
      <td>0.294</td>
      <td>0.195</td>
      <td>0.029</td>
      <td>0.557</td>
      <td>0.528</td>
      <td>30.252172</td>
      <td>16.844959</td>
      <td>37.272727</td>
    </tr>
    <tr>
      <th>theta[10]</th>
      <td>0.378</td>
      <td>0.276</td>
      <td>0.034</td>
      <td>0.703</td>
      <td>0.669</td>
      <td>38.330876</td>
      <td>21.657805</td>
      <td>2.727273</td>
    </tr>
    <tr>
      <th>theta[11]</th>
      <td>0.435</td>
      <td>0.325</td>
      <td>0.014</td>
      <td>0.804</td>
      <td>0.790</td>
      <td>45.263666</td>
      <td>24.923664</td>
      <td>0.909091</td>
    </tr>
    <tr>
      <th>theta[12]</th>
      <td>0.215</td>
      <td>0.142</td>
      <td>0.032</td>
      <td>0.438</td>
      <td>0.406</td>
      <td>23.262086</td>
      <td>12.318593</td>
      <td>66.363636</td>
    </tr>
    <tr>
      <th>theta[13]</th>
      <td>0.224</td>
      <td>0.150</td>
      <td>0.011</td>
      <td>0.421</td>
      <td>0.410</td>
      <td>23.491270</td>
      <td>12.834255</td>
      <td>61.818182</td>
    </tr>
    <tr>
      <th>theta[14]</th>
      <td>0.203</td>
      <td>0.159</td>
      <td>0.006</td>
      <td>0.410</td>
      <td>0.404</td>
      <td>23.147495</td>
      <td>11.631043</td>
      <td>68.181818</td>
    </tr>
    <tr>
      <th>theta[15]</th>
      <td>0.273</td>
      <td>0.207</td>
      <td>0.021</td>
      <td>0.514</td>
      <td>0.493</td>
      <td>28.246819</td>
      <td>15.641748</td>
      <td>39.090909</td>
    </tr>
    <tr>
      <th>theta[16]</th>
      <td>0.210</td>
      <td>0.136</td>
      <td>0.015</td>
      <td>0.380</td>
      <td>0.365</td>
      <td>20.912960</td>
      <td>12.032114</td>
      <td>50.000000</td>
    </tr>
    <tr>
      <th>theta[17]</th>
      <td>0.312</td>
      <td>0.224</td>
      <td>0.020</td>
      <td>0.586</td>
      <td>0.566</td>
      <td>32.429411</td>
      <td>17.876283</td>
      <td>30.909091</td>
    </tr>
    <tr>
      <th>theta[18]</th>
      <td>0.301</td>
      <td>0.227</td>
      <td>0.011</td>
      <td>0.599</td>
      <td>0.588</td>
      <td>33.689918</td>
      <td>17.246030</td>
      <td>28.181818</td>
    </tr>
    <tr>
      <th>theta[19]</th>
      <td>0.280</td>
      <td>0.189</td>
      <td>0.022</td>
      <td>0.530</td>
      <td>0.508</td>
      <td>29.106256</td>
      <td>16.042818</td>
      <td>38.181818</td>
    </tr>
    <tr>
      <th>theta[20]</th>
      <td>0.238</td>
      <td>0.158</td>
      <td>0.031</td>
      <td>0.459</td>
      <td>0.428</td>
      <td>24.522594</td>
      <td>13.636396</td>
      <td>86.363636</td>
    </tr>
    <tr>
      <th>theta[21]</th>
      <td>0.235</td>
      <td>0.157</td>
      <td>0.024</td>
      <td>0.452</td>
      <td>0.428</td>
      <td>24.522594</td>
      <td>13.464508</td>
      <td>73.636364</td>
    </tr>
    <tr>
      <th>theta[22]</th>
      <td>0.235</td>
      <td>0.159</td>
      <td>0.022</td>
      <td>0.455</td>
      <td>0.433</td>
      <td>24.809073</td>
      <td>13.464508</td>
      <td>43.636364</td>
    </tr>
    <tr>
      <th>theta[23]</th>
      <td>0.349</td>
      <td>0.260</td>
      <td>0.006</td>
      <td>0.662</td>
      <td>0.656</td>
      <td>37.586031</td>
      <td>19.996227</td>
      <td>18.181818</td>
    </tr>
    <tr>
      <th>theta[24]</th>
      <td>0.500</td>
      <td>0.398</td>
      <td>0.006</td>
      <td>1.194</td>
      <td>1.188</td>
      <td>68.067386</td>
      <td>28.647890</td>
      <td>13.636364</td>
    </tr>
  </tbody>
</table>
</div>



What about the on-axis source-levels?


```python
onaxis_summary = r.loc['on_axis[0]':'on_axis[24]',:]
onaxis_summary['hdi_width'] = onaxis_summary['hdi_94.5%']-onaxis_summary['hdi_5.5%']
onaxis_summary
```

    <ipython-input-57-bf186442427d>:2: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      onaxis_summary['hdi_width'] = onaxis_summary['hdi_94.5%']-onaxis_summary['hdi_5.5%']





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_5.5%</th>
      <th>hdi_94.5%</th>
      <th>hdi_width</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>on_axis[0]</th>
      <td>99.782</td>
      <td>0.819</td>
      <td>98.475</td>
      <td>101.066</td>
      <td>2.591</td>
    </tr>
    <tr>
      <th>on_axis[1]</th>
      <td>99.908</td>
      <td>0.903</td>
      <td>98.478</td>
      <td>101.368</td>
      <td>2.890</td>
    </tr>
    <tr>
      <th>on_axis[2]</th>
      <td>99.752</td>
      <td>0.863</td>
      <td>98.328</td>
      <td>101.048</td>
      <td>2.720</td>
    </tr>
    <tr>
      <th>on_axis[3]</th>
      <td>99.621</td>
      <td>0.839</td>
      <td>98.174</td>
      <td>100.858</td>
      <td>2.684</td>
    </tr>
    <tr>
      <th>on_axis[4]</th>
      <td>99.998</td>
      <td>0.806</td>
      <td>98.719</td>
      <td>101.279</td>
      <td>2.560</td>
    </tr>
    <tr>
      <th>on_axis[5]</th>
      <td>99.875</td>
      <td>0.855</td>
      <td>98.525</td>
      <td>101.241</td>
      <td>2.716</td>
    </tr>
    <tr>
      <th>on_axis[6]</th>
      <td>99.885</td>
      <td>0.836</td>
      <td>98.638</td>
      <td>101.235</td>
      <td>2.597</td>
    </tr>
    <tr>
      <th>on_axis[7]</th>
      <td>100.046</td>
      <td>0.882</td>
      <td>98.601</td>
      <td>101.315</td>
      <td>2.714</td>
    </tr>
    <tr>
      <th>on_axis[8]</th>
      <td>99.897</td>
      <td>0.856</td>
      <td>98.565</td>
      <td>101.277</td>
      <td>2.712</td>
    </tr>
    <tr>
      <th>on_axis[9]</th>
      <td>99.846</td>
      <td>0.876</td>
      <td>98.416</td>
      <td>101.124</td>
      <td>2.708</td>
    </tr>
    <tr>
      <th>on_axis[10]</th>
      <td>100.044</td>
      <td>0.824</td>
      <td>98.756</td>
      <td>101.386</td>
      <td>2.630</td>
    </tr>
    <tr>
      <th>on_axis[11]</th>
      <td>99.868</td>
      <td>0.851</td>
      <td>98.478</td>
      <td>101.112</td>
      <td>2.634</td>
    </tr>
    <tr>
      <th>on_axis[12]</th>
      <td>100.114</td>
      <td>0.814</td>
      <td>98.935</td>
      <td>101.516</td>
      <td>2.581</td>
    </tr>
    <tr>
      <th>on_axis[13]</th>
      <td>99.792</td>
      <td>0.859</td>
      <td>98.484</td>
      <td>101.204</td>
      <td>2.720</td>
    </tr>
    <tr>
      <th>on_axis[14]</th>
      <td>99.887</td>
      <td>0.859</td>
      <td>98.528</td>
      <td>101.258</td>
      <td>2.730</td>
    </tr>
    <tr>
      <th>on_axis[15]</th>
      <td>100.012</td>
      <td>0.863</td>
      <td>98.609</td>
      <td>101.232</td>
      <td>2.623</td>
    </tr>
    <tr>
      <th>on_axis[16]</th>
      <td>99.592</td>
      <td>0.844</td>
      <td>98.107</td>
      <td>100.814</td>
      <td>2.707</td>
    </tr>
    <tr>
      <th>on_axis[17]</th>
      <td>99.738</td>
      <td>0.812</td>
      <td>98.464</td>
      <td>101.011</td>
      <td>2.547</td>
    </tr>
    <tr>
      <th>on_axis[18]</th>
      <td>99.730</td>
      <td>0.872</td>
      <td>98.309</td>
      <td>101.106</td>
      <td>2.797</td>
    </tr>
    <tr>
      <th>on_axis[19]</th>
      <td>99.837</td>
      <td>0.807</td>
      <td>98.704</td>
      <td>101.180</td>
      <td>2.476</td>
    </tr>
    <tr>
      <th>on_axis[20]</th>
      <td>99.878</td>
      <td>0.869</td>
      <td>98.636</td>
      <td>101.337</td>
      <td>2.701</td>
    </tr>
    <tr>
      <th>on_axis[21]</th>
      <td>99.855</td>
      <td>0.888</td>
      <td>98.413</td>
      <td>101.199</td>
      <td>2.786</td>
    </tr>
    <tr>
      <th>on_axis[22]</th>
      <td>99.894</td>
      <td>0.804</td>
      <td>98.689</td>
      <td>101.257</td>
      <td>2.568</td>
    </tr>
    <tr>
      <th>on_axis[23]</th>
      <td>100.035</td>
      <td>0.826</td>
      <td>98.682</td>
      <td>101.307</td>
      <td>2.625</td>
    </tr>
    <tr>
      <th>on_axis[24]</th>
      <td>99.869</td>
      <td>0.841</td>
      <td>98.468</td>
      <td>101.081</td>
      <td>2.613</td>
    </tr>
  </tbody>
</table>
</div>



### Does the upper 90%ile issue help with source-level estimation?
We've seen that when we have recorded calls from +/- 90$^{\circ}$ of the on-axis direction, the mean SL has been estimated correctly, even though there is still an okay interval of about 3 dB. We can see however, that because any observed value can be reproduced by a very large combination of on-axis source-level and angle. The $\theta$ estimate intervals are *pretty* wide (~$\pm 40^{\circ}$) 

Even if $\theta$ can't be estimated all that well, does choosing the loudest-calls help in decreasing the highest density interval of the on-axis source-level?


```python
louder_sls = measured_SL[measured_SL>np.quantile(measured_SL,0.9)]
```


```python
louder_sls
```




    array([100.16627914,  99.97834302,  98.75851708])



Here, in this case we are kind of lucky as the 90%ile has the true on-axis source level range in it. But we may not always be so lucky, especially if we have a smaller sample size. 


```python
# simulate choosing the loudest calls - only the upper 10%ile of all calls are chosen 
# and the on-axis source-level estimated

m00 = pm.Model()
with m00:
    theta = pm.Uniform('theta',lower=0,upper=np.pi/2,shape=louder_sls.size)
    on_axis = pm.Normal('on_axis',mu=100,sigma=1,shape=louder_sls.size)
    sim = pm.Simulator('sim',fn_measured_level,params=(on_axis,theta),
                                   observed=louder_sls, epsilon=10)
    # -- delete the parallel=True as it always caused issues initially
    trace_lv = pm.sample_smc(draws=500,kernel="ABC")
    idata_lv = az.from_pymc3(trace_lv)
    
```

    Initializing SMC sampler...
    Sampling 2 chains in 2 jobs
    Stage:   0 Beta: 0.439
    Stage:   1 Beta: 1.000
    Stage:   0 Beta: 0.436
    Stage:   1 Beta: 1.000



```python
np.degrees(0.5)
```




    28.64788975654116




```python
loudest_summary = az.summary(idata_lv,kind='stats')
loudest_summary['hdi_width'] = np.abs(loudest_summary['hdi_3%']-loudest_summary['hdi_97%'])
loudest_summary
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_3%</th>
      <th>hdi_97%</th>
      <th>hdi_width</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>on_axis[0]</th>
      <td>99.977</td>
      <td>0.957</td>
      <td>98.227</td>
      <td>101.716</td>
      <td>3.489</td>
    </tr>
    <tr>
      <th>on_axis[1]</th>
      <td>99.964</td>
      <td>1.017</td>
      <td>98.037</td>
      <td>101.827</td>
      <td>3.790</td>
    </tr>
    <tr>
      <th>on_axis[2]</th>
      <td>99.930</td>
      <td>1.025</td>
      <td>97.866</td>
      <td>101.660</td>
      <td>3.794</td>
    </tr>
    <tr>
      <th>theta[0]</th>
      <td>0.460</td>
      <td>0.385</td>
      <td>0.002</td>
      <td>1.359</td>
      <td>1.357</td>
    </tr>
    <tr>
      <th>theta[1]</th>
      <td>0.447</td>
      <td>0.367</td>
      <td>0.002</td>
      <td>1.302</td>
      <td>1.300</td>
    </tr>
    <tr>
      <th>theta[2]</th>
      <td>0.433</td>
      <td>0.370</td>
      <td>0.001</td>
      <td>1.313</td>
      <td>1.312</td>
    </tr>
  </tbody>
</table>
</div>




```python
np.degrees(1.2)
```




    68.75493541569878



### Choosing the loudest calls does not lead to much better estimates of source-level



```python
# simulate the case where all the loudest calls are actually right censored 
# this leaves us with <= 75%ile of all calls. 
fainter_sls = measured_SL[measured_SL<np.quantile(measured_SL,0.75)]

m00b = pm.Model()
with m00b:
    theta = pm.Uniform('theta',lower=0,upper=np.pi/2,shape=fainter_sls.size)
    on_axis = pm.Normal('on_axis',mu=100,sigma=1,shape=fainter_sls.size)
    sim = pm.Simulator('sim',fn_measured_level,params=(on_axis,theta),
                                   observed=fainter_sls, epsilon=10)
    # -- delete the parallel=True as it always caused issues initially
    trace_lvb = pm.sample_smc(draws=500,kernel="ABC")
    idata_lvb = az.from_pymc3(trace_lvb)
    mapb = pm.find_MAP()

```

    Initializing SMC sampler...
    Sampling 2 chains in 2 jobs
    Stage:   0 Beta: 0.051
    Stage:   1 Beta: 0.113
    Stage:   2 Beta: 0.180
    Stage:   3 Beta: 0.251
    Stage:   4 Beta: 0.334
    Stage:   5 Beta: 0.438
    Stage:   6 Beta: 0.559
    Stage:   7 Beta: 0.720
    Stage:   8 Beta: 0.962
    Stage:   9 Beta: 1.000
    Stage:   0 Beta: 0.054
    Stage:   1 Beta: 0.116
    Stage:   2 Beta: 0.187
    Stage:   3 Beta: 0.255
    Stage:   4 Beta: 0.335
    Stage:   5 Beta: 0.430
    Stage:   6 Beta: 0.543
    Stage:   7 Beta: 0.715
    Stage:   8 Beta: 0.959
    Stage:   9 Beta: 1.000




<div>
    <style>
        /* Turns off some styling */
        progress {
            /* gets rid of default border in Firefox and Opera. */
            border: none;
            /* Needs to be in here for Safari polyfill so background images work as expected. */
            background-size: auto;
        }
        .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {
            background: #F44336;
        }
    </style>
  <progress value='1' class='' max='1' style='width:300px; height:20px; vertical-align: middle;'></progress>
  100.00% [1/1 00:00<00:00 logp = -24.669, ||grad|| = 0]
</div>



    



```python
fainter_summary = az.summary(idata_lvb,kind='stats')
fainter_summary['hdi_width'] = np.abs(fainter_summary['hdi_3%']-fainter_summary['hdi_97%'])
fainter_summary
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_3%</th>
      <th>hdi_97%</th>
      <th>hdi_width</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>on_axis[0]</th>
      <td>99.806</td>
      <td>0.933</td>
      <td>98.207</td>
      <td>101.594</td>
      <td>3.387</td>
    </tr>
    <tr>
      <th>on_axis[1]</th>
      <td>99.793</td>
      <td>0.981</td>
      <td>98.086</td>
      <td>101.713</td>
      <td>3.627</td>
    </tr>
    <tr>
      <th>on_axis[2]</th>
      <td>99.782</td>
      <td>0.861</td>
      <td>98.194</td>
      <td>101.419</td>
      <td>3.225</td>
    </tr>
    <tr>
      <th>on_axis[3]</th>
      <td>99.301</td>
      <td>0.886</td>
      <td>97.596</td>
      <td>100.874</td>
      <td>3.278</td>
    </tr>
    <tr>
      <th>on_axis[4]</th>
      <td>99.760</td>
      <td>0.887</td>
      <td>98.231</td>
      <td>101.485</td>
      <td>3.254</td>
    </tr>
    <tr>
      <th>on_axis[5]</th>
      <td>99.741</td>
      <td>0.910</td>
      <td>98.040</td>
      <td>101.399</td>
      <td>3.359</td>
    </tr>
    <tr>
      <th>on_axis[6]</th>
      <td>99.859</td>
      <td>0.919</td>
      <td>98.109</td>
      <td>101.480</td>
      <td>3.371</td>
    </tr>
    <tr>
      <th>on_axis[7]</th>
      <td>99.804</td>
      <td>0.920</td>
      <td>98.061</td>
      <td>101.481</td>
      <td>3.420</td>
    </tr>
    <tr>
      <th>on_axis[8]</th>
      <td>99.848</td>
      <td>0.922</td>
      <td>98.010</td>
      <td>101.389</td>
      <td>3.379</td>
    </tr>
    <tr>
      <th>on_axis[9]</th>
      <td>99.887</td>
      <td>0.960</td>
      <td>98.032</td>
      <td>101.629</td>
      <td>3.597</td>
    </tr>
    <tr>
      <th>on_axis[10]</th>
      <td>99.712</td>
      <td>0.867</td>
      <td>98.084</td>
      <td>101.228</td>
      <td>3.144</td>
    </tr>
    <tr>
      <th>on_axis[11]</th>
      <td>99.753</td>
      <td>0.923</td>
      <td>98.177</td>
      <td>101.551</td>
      <td>3.374</td>
    </tr>
    <tr>
      <th>on_axis[12]</th>
      <td>99.765</td>
      <td>0.878</td>
      <td>97.937</td>
      <td>101.218</td>
      <td>3.281</td>
    </tr>
    <tr>
      <th>on_axis[13]</th>
      <td>99.714</td>
      <td>0.892</td>
      <td>97.840</td>
      <td>101.196</td>
      <td>3.356</td>
    </tr>
    <tr>
      <th>on_axis[14]</th>
      <td>99.895</td>
      <td>0.915</td>
      <td>98.260</td>
      <td>101.641</td>
      <td>3.381</td>
    </tr>
    <tr>
      <th>on_axis[15]</th>
      <td>99.732</td>
      <td>0.938</td>
      <td>97.926</td>
      <td>101.414</td>
      <td>3.488</td>
    </tr>
    <tr>
      <th>on_axis[16]</th>
      <td>99.788</td>
      <td>0.931</td>
      <td>98.037</td>
      <td>101.453</td>
      <td>3.416</td>
    </tr>
    <tr>
      <th>on_axis[17]</th>
      <td>99.964</td>
      <td>0.913</td>
      <td>98.175</td>
      <td>101.561</td>
      <td>3.386</td>
    </tr>
    <tr>
      <th>theta[0]</th>
      <td>0.225</td>
      <td>0.164</td>
      <td>0.010</td>
      <td>0.513</td>
      <td>0.503</td>
    </tr>
    <tr>
      <th>theta[1]</th>
      <td>0.236</td>
      <td>0.185</td>
      <td>0.008</td>
      <td>0.561</td>
      <td>0.553</td>
    </tr>
    <tr>
      <th>theta[2]</th>
      <td>0.224</td>
      <td>0.163</td>
      <td>0.008</td>
      <td>0.512</td>
      <td>0.504</td>
    </tr>
    <tr>
      <th>theta[3]</th>
      <td>0.146</td>
      <td>0.101</td>
      <td>0.012</td>
      <td>0.341</td>
      <td>0.329</td>
    </tr>
    <tr>
      <th>theta[4]</th>
      <td>0.238</td>
      <td>0.159</td>
      <td>0.007</td>
      <td>0.528</td>
      <td>0.521</td>
    </tr>
    <tr>
      <th>theta[5]</th>
      <td>0.232</td>
      <td>0.150</td>
      <td>0.009</td>
      <td>0.500</td>
      <td>0.491</td>
    </tr>
    <tr>
      <th>theta[6]</th>
      <td>0.221</td>
      <td>0.167</td>
      <td>0.009</td>
      <td>0.512</td>
      <td>0.503</td>
    </tr>
    <tr>
      <th>theta[7]</th>
      <td>0.231</td>
      <td>0.167</td>
      <td>0.010</td>
      <td>0.535</td>
      <td>0.525</td>
    </tr>
    <tr>
      <th>theta[8]</th>
      <td>0.229</td>
      <td>0.164</td>
      <td>0.003</td>
      <td>0.521</td>
      <td>0.518</td>
    </tr>
    <tr>
      <th>theta[9]</th>
      <td>0.296</td>
      <td>0.220</td>
      <td>0.004</td>
      <td>0.651</td>
      <td>0.647</td>
    </tr>
    <tr>
      <th>theta[10]</th>
      <td>0.226</td>
      <td>0.161</td>
      <td>0.012</td>
      <td>0.532</td>
      <td>0.520</td>
    </tr>
    <tr>
      <th>theta[11]</th>
      <td>0.230</td>
      <td>0.148</td>
      <td>0.019</td>
      <td>0.505</td>
      <td>0.486</td>
    </tr>
    <tr>
      <th>theta[12]</th>
      <td>0.203</td>
      <td>0.158</td>
      <td>0.007</td>
      <td>0.498</td>
      <td>0.491</td>
    </tr>
    <tr>
      <th>theta[13]</th>
      <td>0.227</td>
      <td>0.184</td>
      <td>0.010</td>
      <td>0.511</td>
      <td>0.501</td>
    </tr>
    <tr>
      <th>theta[14]</th>
      <td>0.225</td>
      <td>0.156</td>
      <td>0.004</td>
      <td>0.494</td>
      <td>0.490</td>
    </tr>
    <tr>
      <th>theta[15]</th>
      <td>0.225</td>
      <td>0.151</td>
      <td>0.012</td>
      <td>0.499</td>
      <td>0.487</td>
    </tr>
    <tr>
      <th>theta[16]</th>
      <td>0.240</td>
      <td>0.159</td>
      <td>0.023</td>
      <td>0.535</td>
      <td>0.512</td>
    </tr>
    <tr>
      <th>theta[17]</th>
      <td>0.255</td>
      <td>0.213</td>
      <td>0.011</td>
      <td>0.568</td>
      <td>0.557</td>
    </tr>
  </tbody>
</table>
</div>




```python
mapb
```




    {'theta_interval__': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
            0.]),
     'on_axis': array([100., 100., 100., 100., 100., 100., 100., 100., 100., 100., 100.,
            100., 100., 100., 100., 100., 100., 100.]),
     'theta': array([0.78539816, 0.78539816, 0.78539816, 0.78539816, 0.78539816,
            0.78539816, 0.78539816, 0.78539816, 0.78539816, 0.78539816,
            0.78539816, 0.78539816, 0.78539816, 0.78539816, 0.78539816,
            0.78539816, 0.78539816, 0.78539816])}




```python
np.degrees(0.785)
```




    44.97718691776962



### Using fainter-calls results in the same interval width as using the loudest calls. However, the $\theta$ estimation interval is much more reduced. 

## Using beam-shapes at different frequencies to estimate on-axis source-levels
The challenge in estimating source-level using the beam-shape for one frequency is that the same measured source-level can be the result of multiple on-axis levels and angle of recording. What if we try to estimate on-axis levels by using two or more beam-shapes? Two or more beam-shapes come about in a call when there is either a harmonic or when the call itself is frequency-modulated (FM). 


```python
multi_rel_levels = np.zeros((theta_values.size,3))
multi_rel_levels[:,0] = theta_values
multi_rel_levels[:,1] = db(piston_in_infinite_baffle(theta_values)) # beam shape for 50kHz
multi_rel_levels[:,2] = db(piston_in_infinite_baffle(theta_values,k=k_value2))# beam shape for 80 kHz

```


```python
plt.figure()
ax0 = plt.subplot(111,projection='polar')
ax0.set_theta_zero_location("N")  # theta=0 at the top
plt.plot(multi_rel_levels[:,0],multi_rel_levels[:,1],label='50kHz')
plt.plot(multi_rel_levels[:,0],multi_rel_levels[:,2],label='80kHz')
plt.legend()
```

The 80kHz beam-shape is more directional, and shows more variation with off-axis angle. Let's see if including this other aspect of information helps us estimate the source-level and off-axis angle better. 

*We are assuming that the sound pressure levels are the same at both frequencies. This may not necessarily be true in case of harmonics, but is likely to be true in case of an FM sweep. Here we are simulating something that matches a typical FM micro-bat call.*


```python
# choose a bunch of random angles of recording and recreate beam-shapes
on_axis_SL = 100 # dB SPL re 20muPa at 1m -   the mean source-level
SL_sd = 1 # a +/- 1dB call-to-call standard deviation

num_points = random_thetas.size
random_rows = np.random.choice(np.arange(multi_rel_levels.shape[0]),num_points)
random_angles = multi_rel_levels[random_rows,0]
off_axis_reduction = multi_rel_levels[random_rows,1:]
on_axis_levels = np.random.normal(on_axis_SL, SL_sd, num_points)
multifreq_measures = off_axis_reduction+np.column_stack((on_axis_levels,on_axis_levels))
```


```python
plt.figure()
plt.plot(multifreq_measures[:,0],label='50khz')
plt.plot(multifreq_measures[:,1],label='80khz')
plt.plot(multifreq_measures[:,0]-multifreq_measures[:,1],label='Difference')
plt.legend();plt.ylabel('Apparent source level (dB SPL re 20$\mu$Pa)');plt.xlabel('Random off-axis angles')

```

We see that there can be very huge differences in the apparent source-level as measured at various angles. For instance, the difference can be upto 40 dB! One lesson I can pick-up from this plot is that perhaps it's better to measure source-levels from lower frequencies than higher frequencies using peak-amplitude etc.

So, there can be strong differences in the apparent source-levels across different frequencies. This might actually help us estimate the $\theta$ and on-axis source-level better?


```python
def fn_measured_level_2freq(on_axis,theta):
    '''Calculates the expected apparent source-level 
    for 50khz and 80kHz given the default piston size.
    '''
    f50khz = on_axis-db(piston_in_infinite_baffle(theta))
    f80khz = on_axis-db(piston_in_infinite_baffle(theta,k=k_value2))
    
    return np.column_stack((f50khz, f80khz))
```


```python
m01 = pm.Model()
with m01:
    theta = pm.Uniform('theta',lower=0,upper=np.pi/2,shape=num_points)
    on_axis = pm.Normal('on_axis',mu=100,sigma=1,shape=num_points)
    sim = pm.Simulator('sim',fn_measured_level_2freq,params=(on_axis,theta),
                                   observed=multifreq_measures, epsilon=10)
    # -- delete the parallel=True as it always caused issues initially
    m01_trace = pm.sample_smc(draws=500,kernel="ABC")
    #idata_lv = az.from_pymc3(trace_lv)
```

    Initializing SMC sampler...
    Sampling 2 chains in 2 jobs
    Stage:   0 Beta: 0.019
    Stage:   1 Beta: 0.039
    Stage:   2 Beta: 0.058
    Stage:   3 Beta: 0.078
    Stage:   4 Beta: 0.098
    Stage:   5 Beta: 0.123
    Stage:   6 Beta: 0.149
    Stage:   7 Beta: 0.179
    Stage:   8 Beta: 0.214
    Stage:   9 Beta: 0.263
    Stage:  10 Beta: 0.319
    Stage:  11 Beta: 0.391
    Stage:  12 Beta: 0.485
    Stage:  13 Beta: 0.614
    Stage:  14 Beta: 0.784
    Stage:  15 Beta: 1.000
    Stage:   0 Beta: 0.018
    Stage:   1 Beta: 0.037
    Stage:   2 Beta: 0.059
    Stage:   3 Beta: 0.079
    Stage:   4 Beta: 0.101
    Stage:   5 Beta: 0.124
    Stage:   6 Beta: 0.149
    Stage:   7 Beta: 0.179
    Stage:   8 Beta: 0.213
    Stage:   9 Beta: 0.256
    Stage:  10 Beta: 0.304
    Stage:  11 Beta: 0.369
    Stage:  12 Beta: 0.451
    Stage:  13 Beta: 0.560
    Stage:  14 Beta: 0.704
    Stage:  15 Beta: 0.881
    Stage:  16 Beta: 1.000



```python
mtrace_summary = az.summary(m01_trace,kind='stats',hdi_prob=0.89)

```

    /home/autumn/anaconda3/envs/statsre/lib/python3.8/site-packages/arviz/data/io_pymc3.py:88: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
      warnings.warn(


### Adding another beam-shape decreases $\theta$ estimation intervals



```python
theta_summary = mtrace_summary.loc['theta[0]':,:]
theta_summary['hdi_width'] = theta_summary['hdi_94.5%']-theta_summary['hdi_5.5%']
theta_summary['hdiwidth_deg'] = np.degrees(theta_summary['hdi_width'])
theta_summary
```

    <ipython-input-31-245cfa5cc629>:2: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      theta_summary['hdi_width'] = theta_summary['hdi_94.5%']-theta_summary['hdi_5.5%']
    <ipython-input-31-245cfa5cc629>:3: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      theta_summary['hdiwidth_deg'] = np.degrees(theta_summary['hdi_width'])





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_5.5%</th>
      <th>hdi_94.5%</th>
      <th>hdi_width</th>
      <th>hdiwidth_deg</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>theta[0]</th>
      <td>0.125</td>
      <td>0.071</td>
      <td>0.019</td>
      <td>0.223</td>
      <td>0.204</td>
      <td>11.688339</td>
    </tr>
    <tr>
      <th>theta[1]</th>
      <td>0.111</td>
      <td>0.075</td>
      <td>0.009</td>
      <td>0.207</td>
      <td>0.198</td>
      <td>11.344564</td>
    </tr>
    <tr>
      <th>theta[2]</th>
      <td>0.118</td>
      <td>0.064</td>
      <td>0.025</td>
      <td>0.204</td>
      <td>0.179</td>
      <td>10.255945</td>
    </tr>
    <tr>
      <th>theta[3]</th>
      <td>0.131</td>
      <td>0.067</td>
      <td>0.031</td>
      <td>0.223</td>
      <td>0.192</td>
      <td>11.000790</td>
    </tr>
    <tr>
      <th>theta[4]</th>
      <td>0.171</td>
      <td>0.124</td>
      <td>0.021</td>
      <td>0.342</td>
      <td>0.321</td>
      <td>18.391945</td>
    </tr>
    <tr>
      <th>theta[5]</th>
      <td>0.216</td>
      <td>0.151</td>
      <td>0.018</td>
      <td>0.402</td>
      <td>0.384</td>
      <td>22.001579</td>
    </tr>
    <tr>
      <th>theta[6]</th>
      <td>0.226</td>
      <td>0.161</td>
      <td>0.026</td>
      <td>0.422</td>
      <td>0.396</td>
      <td>22.689129</td>
    </tr>
    <tr>
      <th>theta[7]</th>
      <td>0.134</td>
      <td>0.074</td>
      <td>0.020</td>
      <td>0.236</td>
      <td>0.216</td>
      <td>12.375888</td>
    </tr>
    <tr>
      <th>theta[8]</th>
      <td>0.137</td>
      <td>0.090</td>
      <td>0.012</td>
      <td>0.260</td>
      <td>0.248</td>
      <td>14.209353</td>
    </tr>
    <tr>
      <th>theta[9]</th>
      <td>0.130</td>
      <td>0.087</td>
      <td>0.016</td>
      <td>0.251</td>
      <td>0.235</td>
      <td>13.464508</td>
    </tr>
    <tr>
      <th>theta[10]</th>
      <td>0.123</td>
      <td>0.067</td>
      <td>0.020</td>
      <td>0.211</td>
      <td>0.191</td>
      <td>10.943494</td>
    </tr>
    <tr>
      <th>theta[11]</th>
      <td>0.138</td>
      <td>0.073</td>
      <td>0.025</td>
      <td>0.241</td>
      <td>0.216</td>
      <td>12.375888</td>
    </tr>
    <tr>
      <th>theta[12]</th>
      <td>0.118</td>
      <td>0.067</td>
      <td>0.025</td>
      <td>0.208</td>
      <td>0.183</td>
      <td>10.485128</td>
    </tr>
    <tr>
      <th>theta[13]</th>
      <td>0.147</td>
      <td>0.079</td>
      <td>0.026</td>
      <td>0.254</td>
      <td>0.228</td>
      <td>13.063438</td>
    </tr>
    <tr>
      <th>theta[14]</th>
      <td>0.186</td>
      <td>0.119</td>
      <td>0.028</td>
      <td>0.338</td>
      <td>0.310</td>
      <td>17.761692</td>
    </tr>
    <tr>
      <th>theta[15]</th>
      <td>0.122</td>
      <td>0.072</td>
      <td>0.023</td>
      <td>0.226</td>
      <td>0.203</td>
      <td>11.631043</td>
    </tr>
    <tr>
      <th>theta[16]</th>
      <td>0.144</td>
      <td>0.088</td>
      <td>0.020</td>
      <td>0.258</td>
      <td>0.238</td>
      <td>13.636396</td>
    </tr>
    <tr>
      <th>theta[17]</th>
      <td>0.124</td>
      <td>0.068</td>
      <td>0.026</td>
      <td>0.220</td>
      <td>0.194</td>
      <td>11.115381</td>
    </tr>
    <tr>
      <th>theta[18]</th>
      <td>0.136</td>
      <td>0.073</td>
      <td>0.035</td>
      <td>0.243</td>
      <td>0.208</td>
      <td>11.917522</td>
    </tr>
    <tr>
      <th>theta[19]</th>
      <td>0.101</td>
      <td>0.059</td>
      <td>0.021</td>
      <td>0.180</td>
      <td>0.159</td>
      <td>9.110029</td>
    </tr>
    <tr>
      <th>theta[20]</th>
      <td>0.118</td>
      <td>0.069</td>
      <td>0.023</td>
      <td>0.215</td>
      <td>0.192</td>
      <td>11.000790</td>
    </tr>
    <tr>
      <th>theta[21]</th>
      <td>0.121</td>
      <td>0.075</td>
      <td>0.019</td>
      <td>0.227</td>
      <td>0.208</td>
      <td>11.917522</td>
    </tr>
    <tr>
      <th>theta[22]</th>
      <td>0.103</td>
      <td>0.068</td>
      <td>0.010</td>
      <td>0.195</td>
      <td>0.185</td>
      <td>10.599719</td>
    </tr>
    <tr>
      <th>theta[23]</th>
      <td>0.127</td>
      <td>0.075</td>
      <td>0.011</td>
      <td>0.228</td>
      <td>0.217</td>
      <td>12.433184</td>
    </tr>
    <tr>
      <th>theta[24]</th>
      <td>0.115</td>
      <td>0.074</td>
      <td>0.014</td>
      <td>0.217</td>
      <td>0.203</td>
      <td>11.631043</td>
    </tr>
  </tbody>
</table>
</div>



We see that the 89% highest density interval for all $\theta$ values is now reduced to at most 20$^{\circ}$ from a max of 40$^{\circ}$ before with only one beam-shape included. This is quite the drop in uncertainty. 


```python
onaxis_summary = mtrace_summary.loc['on_axis[0]':'on_axis[24]',:]
onaxis_summary['hdi_width'] = onaxis_summary['hdi_94.5%']-onaxis_summary['hdi_5.5%']
onaxis_summary
```

    <ipython-input-32-0addd4271c23>:2: SettingWithCopyWarning: 
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead
    
    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      onaxis_summary['hdi_width'] = onaxis_summary['hdi_94.5%']-onaxis_summary['hdi_5.5%']





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>hdi_5.5%</th>
      <th>hdi_94.5%</th>
      <th>hdi_width</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>on_axis[0]</th>
      <td>99.634</td>
      <td>0.796</td>
      <td>98.397</td>
      <td>100.846</td>
      <td>2.449</td>
    </tr>
    <tr>
      <th>on_axis[1]</th>
      <td>99.555</td>
      <td>0.855</td>
      <td>98.203</td>
      <td>100.848</td>
      <td>2.645</td>
    </tr>
    <tr>
      <th>on_axis[2]</th>
      <td>99.584</td>
      <td>0.881</td>
      <td>98.033</td>
      <td>100.773</td>
      <td>2.740</td>
    </tr>
    <tr>
      <th>on_axis[3]</th>
      <td>99.613</td>
      <td>0.803</td>
      <td>98.290</td>
      <td>100.817</td>
      <td>2.527</td>
    </tr>
    <tr>
      <th>on_axis[4]</th>
      <td>99.763</td>
      <td>0.919</td>
      <td>98.423</td>
      <td>101.290</td>
      <td>2.867</td>
    </tr>
    <tr>
      <th>on_axis[5]</th>
      <td>99.914</td>
      <td>0.873</td>
      <td>98.468</td>
      <td>101.299</td>
      <td>2.831</td>
    </tr>
    <tr>
      <th>on_axis[6]</th>
      <td>99.796</td>
      <td>0.859</td>
      <td>98.421</td>
      <td>101.078</td>
      <td>2.657</td>
    </tr>
    <tr>
      <th>on_axis[7]</th>
      <td>99.820</td>
      <td>0.777</td>
      <td>98.561</td>
      <td>101.055</td>
      <td>2.494</td>
    </tr>
    <tr>
      <th>on_axis[8]</th>
      <td>99.769</td>
      <td>0.842</td>
      <td>98.504</td>
      <td>101.221</td>
      <td>2.717</td>
    </tr>
    <tr>
      <th>on_axis[9]</th>
      <td>99.521</td>
      <td>0.822</td>
      <td>98.263</td>
      <td>100.815</td>
      <td>2.552</td>
    </tr>
    <tr>
      <th>on_axis[10]</th>
      <td>99.746</td>
      <td>0.853</td>
      <td>98.444</td>
      <td>101.135</td>
      <td>2.691</td>
    </tr>
    <tr>
      <th>on_axis[11]</th>
      <td>99.922</td>
      <td>0.803</td>
      <td>98.650</td>
      <td>101.170</td>
      <td>2.520</td>
    </tr>
    <tr>
      <th>on_axis[12]</th>
      <td>99.455</td>
      <td>0.831</td>
      <td>98.067</td>
      <td>100.740</td>
      <td>2.673</td>
    </tr>
    <tr>
      <th>on_axis[13]</th>
      <td>99.497</td>
      <td>0.810</td>
      <td>98.238</td>
      <td>100.807</td>
      <td>2.569</td>
    </tr>
    <tr>
      <th>on_axis[14]</th>
      <td>100.029</td>
      <td>0.796</td>
      <td>98.567</td>
      <td>101.159</td>
      <td>2.592</td>
    </tr>
    <tr>
      <th>on_axis[15]</th>
      <td>99.562</td>
      <td>0.874</td>
      <td>98.182</td>
      <td>100.990</td>
      <td>2.808</td>
    </tr>
    <tr>
      <th>on_axis[16]</th>
      <td>99.805</td>
      <td>0.761</td>
      <td>98.635</td>
      <td>101.051</td>
      <td>2.416</td>
    </tr>
    <tr>
      <th>on_axis[17]</th>
      <td>99.563</td>
      <td>0.866</td>
      <td>98.182</td>
      <td>100.871</td>
      <td>2.689</td>
    </tr>
    <tr>
      <th>on_axis[18]</th>
      <td>99.878</td>
      <td>0.796</td>
      <td>98.725</td>
      <td>101.176</td>
      <td>2.451</td>
    </tr>
    <tr>
      <th>on_axis[19]</th>
      <td>99.321</td>
      <td>0.849</td>
      <td>98.041</td>
      <td>100.722</td>
      <td>2.681</td>
    </tr>
    <tr>
      <th>on_axis[20]</th>
      <td>99.691</td>
      <td>0.819</td>
      <td>98.410</td>
      <td>100.990</td>
      <td>2.580</td>
    </tr>
    <tr>
      <th>on_axis[21]</th>
      <td>99.589</td>
      <td>0.844</td>
      <td>98.252</td>
      <td>100.862</td>
      <td>2.610</td>
    </tr>
    <tr>
      <th>on_axis[22]</th>
      <td>99.143</td>
      <td>0.903</td>
      <td>97.785</td>
      <td>100.678</td>
      <td>2.893</td>
    </tr>
    <tr>
      <th>on_axis[23]</th>
      <td>99.776</td>
      <td>0.808</td>
      <td>98.426</td>
      <td>101.017</td>
      <td>2.591</td>
    </tr>
    <tr>
      <th>on_axis[24]</th>
      <td>99.746</td>
      <td>0.845</td>
      <td>98.352</td>
      <td>101.009</td>
      <td>2.657</td>
    </tr>
  </tbody>
</table>
</div>



There is not much drop in the on-axis level uncertainty. The same ~3 dB interval remains as in the single beam-shape example.  Remember, we have maintained everything the same (sample size, on-axis level$\pm$sd), and the only new thing is the additional beam-shape information. 

## Summary
Measuring on-axis source-levels for bats flying in the field is no easy task. Even if the bat is flying towards the microphones, the call may not be directed at the mics straight away. In this notebook I have shown how explicitly modelling beam-shape effects increases the accuracy of on-axis source level estimation, when we include all calls - and not only the 'loudest' calls. Till now, researchers have used an upper percentile threshold to estimate the on-axis source levels. Using the approach outlined in this notebook, we can thus actually estimate the on-axis levels even with relatively faint calls.  


```python
%load_ext watermark
%watermark -n -u -v -iv `b -w
```

    Last updated: Sun Mar 14 2021
    
    Python implementation: CPython
    Python version       : 3.8.5
    IPython version      : 7.19.0
    
    numpy     : 1.19.5
    arviz     : 0.11.0
    pymc3     : 3.11.0
    matplotlib: 3.3.4
    pandas    : 1.2.1
    
    Watermark: 2.1.0
    



```python
import os

```


```python
os.getcwd()
```




    '/home/autumn/Documents/trying_out/me-statistical-rethinking/paramsims'




```python

```

a8bf2461f2ffbe3d16b5f56e4c338e36d6c87323
ipynb
Jupyter Notebook
paramsims/source-level_estim.ipynb
thejasvibr/tbr-statisticalrethinking
cffda510ec3d7d15c603a89482b1e5cb4f7a3a4a
MIT
paramsims/source-level_estim.ipynb
thejasvibr/tbr-statisticalrethinking
cffda510ec3d7d15c603a89482b1e5cb4f7a3a4a
MIT
paramsims/source-level_estim.ipynb
thejasvibr/tbr-statisticalrethinking
cffda510ec3d7d15c603a89482b1e5cb4f7a3a4a
MIT
Qwen/Qwen-72B
1. YES
2. YES