fits the KDE to the data by estimating the bandwidths and storing the data
			
			Parameters
			----------
				data: 2d-array, shape N x M
					N datapoints in an M dimensional space to which the KDE is fit
				weights: 1d array
					N weights, one for every data point.
					They will be normalized to sum up to one
				fix_boundary_effects: bool
					whether to reweigh points close to the bondary no fix the pdf
				bw_estimator: str
					allowed values are 'scott' and 'mlcv' for Scott's rule of thumb
					and the maximum likelihood via cross-validation
				efficient_bw_estimation: bool
					if true, start bandwidth optimization from the previous value, otherwise
					start from Scott's values
				update_bandwidths: bool
					whether to update the bandwidths at all

def fit(self, data, weights=None, bw_estimator='scott', efficient_bw_estimation=True, update_bandwidth=True):
		"""
			fits the KDE to the data by estimating the bandwidths and storing the data
			
			Parameters
			----------
				data: 2d-array, shape N x M
					N datapoints in an M dimensional space to which the KDE is fit
				weights: 1d array
					N weights, one for every data point.
					They will be normalized to sum up to one
				fix_boundary_effects: bool
					whether to reweigh points close to the bondary no fix the pdf
				bw_estimator: str
					allowed values are 'scott' and 'mlcv' for Scott's rule of thumb
					and the maximum likelihood via cross-validation
				efficient_bw_estimation: bool
					if true, start bandwidth optimization from the previous value, otherwise
					start from Scott's values
				update_bandwidths: bool
					whether to update the bandwidths at all
		"""
		
		if self.data is None:
			# overwrite some values in case this is the first fit of the KDE
			efficient_bw_estimation = False
			update_bandwidth=True

		self.data = np.asfortranarray(data)
		for i,k in enumerate(self.kernels):
				self.kernels[i].data = self.data[:,i]

		self.weights = self._normalize_weights(weights)

		if not update_bandwidth:
			return
		
		if not efficient_bw_estimation or bw_estimator == 'scott':
			# inspired by the the statsmodels code
			sigmas = np.std(self.data, ddof=1, axis=0)
			IQRs = np.subtract.reduce(np.percentile(self.data, [75,25], axis=0))
			self.bandwidths = 1.059 * np.minimum(sigmas, IQRs) * np.power(self.data.shape[0], -0.2)
			# crop bandwidths for categorical parameters
			self.bandwidths = np.clip(self.bandwidths , self.min_bandwidth, self.bw_clip)
			
		if bw_estimator == 'mlcv':
			# optimize bandwidths here
			def opt_me(bw):
				self.bandwidths=bw
				self._set_kernel_bandwidths()
				return(self.loo_negloglikelihood())
			
			res = spo.minimize(opt_me, self.bandwidths, bounds=self.bw_bounds, method='SLSQP')
			self.optimizer_result = res
			self.bandwidths[:] = res.x
		self._set_kernel_bandwidths()