where \\(m(\\vec{x})\\) is the mean function and \\(K(\\vec{x},\\vec{x})\\) represents the covariance between function values. For this work, a regression mean function is assumed. Using this form, the mean function has the following form:
\n where \\(h^{T}(\\vec{x})\\) represents a column vector containing the basis functions of the basis evaluated at the points \\(\\vec{x}\\). The regression parameters \\(\\beta\\) are treated as part of the Kriging model and are determined while constructioning the model. Using this form of the mean function yields a Universal Kriging model. The case of a (polyorder=0) regression (where the vector \\(h(\\vec{x})\\) reduces to unity) is referred to as Ordinary Kriging and is also covered by this functional form. The assumption of a vague prior on the regression parameters gives the following closed form for the parameters:
\n where \\(K\\) is the covariance matrix between the training data. For a Kriging model, the covariance between function values is assumed to be only a function of the distance between points. The multi-dimension covariance function is constructed using a tensor product of one dimension functions. The multi-dimension covariance is calculated in fit(), which calls the static method covfn. The elements in ths covariance matrix are given as:
\n\n\\[ K_{i,j} = cov(y_{i},y_{j}) = \\sigma^{2} k(\\vec{X}_{i},\\vec{X}_{j}; \\theta) + \\sigma^{2}_{n} \\delta_{i,j} \\]\n
\n The parameters \\(\\sigma\\) and \\(\\theta\\) (and in some cases \\(\\sigma_{n}\\)) are denoted as hyperparameters and are determined maximizing the likelihood equation for the Kriging model. This likelihood gives the probability that a Gaussian process with specified hyperparamters describes the training data X and Y. By picking the hyperparameters that maximize this probability, a Kriging model that best describes the data can be constructed. The hyperparameters are determined based on the likelihood equation for a gaussian process with a vague prior on the regression parameters. This likelihood is computed in function llhobj which uses the Scipy Module, minimize. Using this equation, optimization is used to determine all of the parameters, including the covariance magnitude \\(\\sigma\\) and noise \\(\\sigma_{n}\\). This way of determining hyperparameters should be used when the noise level of the function needs to be fitted.
\n With the regression and covariance parameters determined, the final processed data can be constructed using the inverse of the covariance matrix. To make predictions from the Kriging model, the following vector is required:
\n where \\( K \\) is the covariance matrix, the product \\(H^{T} \\beta\\) represents the mean function evaluated at the training points and \\(Y\\) represents the function values at the training points. Using this processed data, the regression parameters and covariance parameters & predictions can be made.\n\n
\n \n
| | Model predictions throughout the domain are determined by sampling from the conditional distribution \\(y_* | \\vec{X},Y\\) using the covariance between points in the domain where \\(\\vec{X},Y\\) are the input and output training data. The posterior mean predictions for an explicit mean are given by the formula: \n \n\\[ y(\\vec{x}_{*}) | \\vec{X},Y,m(x) = m(\\vec{x}_{*}) + k_*^T K^{-1} (Y-m(\\vec{x}_{*})) \\]\n \n where \\(k_{*}^{T}\\) represents the covariance between the test point, \\(\\vec{x}_{*}\\), and the training points \\(\\vec{X}\\) (a row vector of length ntot).
\n For a regression mean function, the function predictions take the form of: \n\n\\[ y(\\vec{x}_{*}) | \\vec{X},Y,\\beta = h^{T}(\\vec{x}_{*}) \\beta + k_*^T K^{-1} (Y-H^{T} \\beta) \\]\n \n The regression parameters \\(\\beta\\) and the hyperparamters in the covariance function are supplied by fit(). Using only this data, function predictions can be made; however, the construction and inverse of the covariance matrix can make the function predictions expensive. Because this matrix is inverted during the construction of the Kriging model, this work can be re-used for function predictions. Defining the processed data \\(V\\) as: \n\n\\[ V = K^{-1} (Y - H^{T} \\beta) \\]\n \n the function predictions are given by: \n\n\\[ y(\\vec{x}_{*}) | \\vec{X},Y,\\beta = h^{T}(\\vec{x}_{*}) \\beta + k_*^T V \\]\n ","metadata":{}},{"cell_type":"code","source":"from sklearn.base import BaseEstimator,RegressorMixin\nfrom numpy.linalg import cholesky, det, lstsq, inv, pinv\nfrom scipy.optimize import minimize\nimport numpy as np\nfrom sklearn.preprocessing import PolynomialFeatures\npi = 4.0*np.arctan(1.0)\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\n# Universal Kriging Model (Polynomial Regression + Full Gaussian Process Regression Model)\n# Commonly used Ensemble Approach for geospatial interpolation and \n\nclass Kriging(BaseEstimator,RegressorMixin):\n \n def __init__(self,kernel='rbf',theta=10.0,sigma=10.0,sigma_n=1,opt=True,polyorder=2):\n self.theta = theta\n self.sigma = sigma\n self.sigma_n = sigma_n\n self.opt = opt\n self.polyorder = polyorder \n Kriging.kernel = kernel \n\n ''' local covariance functions '''\n @staticmethod\n def covfn(X0,X1,theta=1.0,sigma=1.0):\n\n ''' Radial Basis Covariance Function '''\n if(Kriging.kernel == 'rbf'):\n r = np.sum(X0**2,1).reshape(-1,1) + np.sum(X1**2,1) - 2 * np.dot(X0,X1.T)\n return sigma**2 * np.exp(-0.5/theta**2*r)\n\n ''' Matern Covariance Class of Funtions '''\n if(Kriging.kernel == 'matern'):\n lid=1\n r = np.sum(X0**2,1)[:,None] + np.sum(X1**2,1) - 2 * np.dot(X0,X1.T)\n if(lid==1):\n return sigma**2 * np.exp(-r/theta)\n elif(lid==2):\n ratio = r/theta\n v1 = (1.0+np.sqrt(3)*ratio)\n v2 = np.exp(-np.sqrt(3)*ratio)\n return sigma**2*v1*v2\n elif(lid==3):\n ratio = r/theta\n v1 = (1.0+np.sqrt(5)*ratio+(5.0/3.0)*ratio**2)\n v2 = np.exp(-np.sqrt(5)*ratio)\n return sigma**2*v1*v2\n else:\n print('Covariance Function not defined')\n \n ''' Train the Model'''\n def fit(self,X,y):\n \n ''' Working w/ numpy matrices'''\n if(type(X) is np.ndarray):\n self.X = X;self.y = y\n else:\n self.X = X.values; self.y = y.values\n self.ntot,ndim = self.X.shape\n \n # Collocation Matrix\n self.poly = PolynomialFeatures(self.polyorder)\n self.H = self.poly.fit_transform(self.X)\n \n ''' Optimisation Objective Function '''\n # Optimisation of hyperparameters via the objective funciton\n def llhobj(X,y,noise):\n \n # Simplified Variant\n def llh_dir(hypers):\n K = self.covfn(X,X,theta=hypers[0],sigma=hypers[1]) + noise**2 * np.eye(self.ntot)\n return 0.5 * np.log(det(K)) + \\\n 0.5 * y.T.dot(inv(K).dot(y)).ravel()[0] + 0.5 * self.ntot * np.log(2*pi)\n\n # Full Likelihood Equation\n def nll_full(hypers):\n K = self.covfn(X,X,theta=hypers[0],sigma=hypers[1]) + noise**2 * np.eye(self.ntot)\n L = cholesky(K)\n return np.sum(np.log(np.diagonal(L))) + \\\n 0.5 * y.T.dot(lstsq(L.T, lstsq(L,y)[0])[0]) + \\\n 0.5 * self.ntot * np.log(2*pi)\n \n return llh_dir # return one of the two, simplified variant doesn't always work well\n \n ''' Update hyperparameters based on set objective function '''\n if(self.opt==True):\n # define the objective funciton\n objfn = llhobj(self.X,self.y,self.sigma_n)\n # search for the optimal hyperparameters based on given relation\n res = minimize(fun=objfn,x0=[1,1],\n method='Nelder-Mead',tol=1e-6)\n self.theta,self.sigma = res.x # update the hyperparameters to \n\n self.HT = self.H.T\n self.Kmat = self.covfn(self.X,self.X,self.theta,self.sigma) \\\n + self.sigma_n**2 * np.eye(self.ntot) # Covariance Matrix (Train/Train)\n self.IKmat = pinv(self.Kmat) # Pseudo Matrix Inversion (More Stable)\n\n self.HK = np.dot(self.HT,self.IKmat) # HK^-1\n HKH = np.dot(self.HK,self.H) # HK^-1HT\n self.A = inv(HKH) # Variance-Covariance Weighted LS Matrix\n\n self.W = np.dot(self.IKmat,self.y)\n Q = np.dot(self.HT,self.W)\n self.beta = np.dot(self.A,Q) # Regression coefficients\n self.V = self.W - np.dot(self.IKmat,self.H).dot(self.beta) # K^{-1} (Y - H^{T} * beta)\n \n return self # return class & use w/ predict()\n\n ''' Posterior Prediction; '''\n # Make a prediction based on what the model has learned \n def predict(self,Xm):\n \n ''' Working w/ numpy matrices'''\n if(type(Xm) is np.ndarray):\n self.Xm = Xm\n else:\n self.Xm = Xm.values\n self.mtot,ndim = self.Xm.shape\n \n self.Hm = self.poly.fit_transform(self.Xm) # Collocation Matrix\n self.Kmat = self.covfn(self.X,self.Xm,self.theta,self.sigma) # Covariance Matrix (Train/Test)\n yreg = np.dot(self.Hm,self.beta) # Mean Prediction based on Regression\n ykr = np.dot(self.Kmat.T,self.V) # posterior mean predictions for an explicit mean \n\n return yreg + ykr","metadata":{"execution":{"iopub.status.busy":"2022-02-17T16:12:08.011877Z","iopub.execute_input":"2022-02-17T16:12:08.012168Z","iopub.status.idle":"2022-02-17T16:12:08.045201Z","shell.execute_reply.started":"2022-02-17T16:12:08.01214Z","shell.execute_reply":"2022-02-17T16:12:08.043684Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5 | CREATING SPECTOGRAMS\n\n\n\n- get_spectograms is an extended function of the one posted another notebook **[[Keras, Inference] BirdCLEF2021 starter](https://www.kaggle.com/shtrausslearning/keras-inference-birdclef2021-starter)**\n- The function is used for creating a spectogram for different audio chunks of the split input audio signal.\n- I've included several things in the function that can be useful to visualise during general EDA & we will look through one recording that is split into several segments (all contents for one case only.\n\n**Function Arguments**\n\n- **filepath** : string containg path to audio
\n- **primary_label** : desired output name
\n- **output_dir** : directory in which the spectograms are saved
\n- **save_id** : Save spectogram output or not
\n- **audio_id** : display segment audio recording
\n- **select_id** : training / soundscape data is looked at ( soundscape just adds label )
\n- **plot_id** : output general plots ( those displayed in the next example )\n- **store_id** : save peak values to global loop ( when actually looping through all audios )\n- **height** : plotly figure height","metadata":{}},{"cell_type":"code","source":"def get_spectrograms(filepath=None,primary_label=None, output_dir=None,save_id=False,\n audio_id=False,select_id='train',plot_id=False,store_id=False,height=400):\n\n # duration is set from global variable\n sig, rate = librosa.load(filepath, sr=cfg.sr, offset=cfg.offset, duration=cfg.cutoff) # read audio data\n sig_splits = split_signal(sig) # split the signal into parts \n fig = make_subplots(rows=3, cols=len(sig_splits))\n \n # Extract mel spectrograms for each audio chunk\n s_cnt = 0; path_all = []; jj=0; kk= cfg.offset; lst_peak_model = []\n for chunk in sig_splits:\n \n kk=+ cfg.offset + cfg.sl\n \n # Play Audio\n if(audio_id):\n print(f'Audio ID: {kk}')\n display(ipd.display(ipd.Audio(data=chunk, rate=cfg.sr)))\n \n kk+=5; jj+=1;\n mel = librosa.feature.melspectrogram(chunk,sr=cfg.sr,\n fmin=cfg.fmin, \n fmax=cfg.fmax,\n n_mels=cfg.sshape[0],\n n_fft = cfg.nfft) \n mel_spec = librosa.power_to_db(mel**cfg.pw, ref=np.min)\n mel_spec_disp = mel_spec.copy()\n mel_spec_disp -= mel_spec_disp.min(); mel_spec_disp /= mel_spec_disp.max() # numpy format\n if(plot_id):\n fig.add_trace(go.Heatmap(z=mel_spec,colorscale='viridis',showscale=False),1,jj)\n \n # add sandscape results annotation\n if(select_id is 'soundscape'):\n data = get_subset()\n record_id = filepath.split('_')[1].split('/')[1]\n get_birds = data.pd_scape[data.pd_scape['audio_id']==int(record_id)][['seconds','birds']]\n get_birds_v= get_birds[get_birds['seconds'] == kk]['birds'].values[0]\n \n '''Get Mel Max/Min/Mean Values (TIME)'''\n ldf = pd.DataFrame(mel_spec)\n maxst = ldf.describe().loc['max',:]\n if(plot_id):\n fig.add_trace(go.Scatter(y=maxst,line=dict(color=lst_col[0])),2,jj) # chunk signal data\n\n '''Get Mel Max/Min/Mean Values (FREQUENCY)'''\n ldf = pd.DataFrame(mel_spec.T)\n maxsf = ldf.describe().loc['max',:]\n if(plot_id):\n fig.add_trace(go.Scatter(y=maxsf,line = dict(color=lst_col[1]),name='spect-max'),3,jj)\n \n ''' TIME BIN MODEL PREDICTION '''\n # Split the data into time bins, splitting data into bins -> broadband model attempt\n\n lst_df = split_dataframe(ldf,chunk_size=cfg.model_bins) \n lst_ensemble = []; ym = 0; ii=0\n for tdf in lst_df:\n ii+=1; maxs = tdf.describe().loc['max',:] # maximum value in frequency band (plot)\n model = Kriging(opt=False,theta=4)\n model.fit(X=tdf.T.index[:,None],y=maxs)\n Xm = np.arange(0,maxsf.shape[0],1)[:,None]\n ym+= model.predict(Xm)\n y_ens = ym/float(ii) # ensemble model\n if(plot_id):\n fig.add_trace(go.Scatter(y=y_ens,line = dict(color=lst_col[0]),name='spect-bin-model'),3,jj)\n \n ''' Find Peaks in Data '''\n # find peaks in bin model \n peaks, _ = find_peaks(y_ens, height=0) # find peaks in ensemble model\n if(store_id):\n glst_peak_model.extend(peaks) \n if(plot_id):\n fig.add_trace(go.Scatter(x=peaks,y=y_ens[peaks],mode='markers',marker=dict(color='black')),3,jj)\n \n ''' (SAVE) MELSPECTOGRAM '''\n if(save_id):\n mel_spec -= mel_spec.min(); mel_spec /= mel_spec.max() # numpy format\n save_dir = os.path.join(output_dir, primary_label)\n if not os.path.exists(save_dir): os.makedirs(save_dir)\n fname = filepath.rsplit(os.sep, 1)[-1].rsplit('.', 1)[0] + '_' + str(s_cnt) + '.png'\n save_path = os.path.join(save_dir,fname)\n im = Image.fromarray(mel_spec * 255.0).convert(\"L\")\n im.save(save_path)\n \n # add filepath to list\n path_all.append(save_path)\n s_cnt += 1\n \n if(plot_id):\n fig.update_layout(margin=dict(l=0, r=0, t=30, b=0),height=height,showlegend=False)\n fig.update_layout(template='plotly_white',font=dict(family='sans-serif',size=14))\n fig.show()\n \n# if(plot_id):\n# fig.update_layout(margin=dict(l=0, r=0, t=30, b=0),height=height,coloraxis_showscale=False,showlegend=False)\n# fig.update_layout(template='plotly_white',font=dict(family='sans-serif',size=14))\n# fig.show()\n \n return path_all # return list of pathways to created spectograms","metadata":{"execution":{"iopub.status.busy":"2022-02-17T16:12:08.047178Z","iopub.execute_input":"2022-02-17T16:12:08.047654Z","iopub.status.idle":"2022-02-17T16:12:08.077008Z","shell.execute_reply.started":"2022-02-17T16:12:08.047603Z","shell.execute_reply":"2022-02-17T16:12:08.075804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n 5.2 | SAMPLE VISUALISATION CASE \n \n 5.4 | SOME OVERSERVATIONS \n \n 5.5 | MODELS OR FASTER METHODS \n \n 5.6 | MISSING POTENTIAL BIRD CALLS \n |