Probabilistic Outlier Detection (edited + clarified)
In essence anomaly detection is about finding a metric with which to measure the similarity between instances and then determining a threshold when crossed constitutes an anomaly. There are parametric anomaly detection algorithms which require some hyper-parameters to be set or information about the distribution to be known. There are also non-parametric techniques which do not have these stringent requirements. 

Let's consider the simple case of identifying anomalies in a population based on height. We can assume this to be a 1D Gaussian distribution. Now we need a similarity metric, let's use Euclidean distance for simplicity (the difference in individual's heights). Next we need a threshold. 

Classical statistics

Classical methods would parametrize the Gaussian distribution of your population. Let us assume that after polling our population we identify that the mean height is 170cm with a variance of 15cm. Now we can set the threshold to be any multiple of the variance $\sigma$. Typically we choose $3\sigma$ as a thresholds since the ingroup would include 99.7% of the population, all those outside of the $3\sigma$ range can thus confidently be said to be anomalous. 

You can also look at the application of a generalized likelihood ratio test (GLRT) to find anomalies.

Let's visualize this

import numpy as np

data = np.random.normal(170, 15, 10000)

plt.hist(data)
plt.plot([215,215],[0,3000])
plt.plot([125,125],[0,3000])
plt.show()




Parametric techniques

The examples will require that we have prior information regarding the distribution of the data. This technique will fit an ellipse around the Gaussian distribution of our data in order to contain the elements with the highest in-group correlation.

from sklearn.covariance import EllipticEnvelope
estimator = EllipticEnvelope()
estimator.fit(data.reshape(-1,1))

plt.scatter(data, np.random.randint(0,1,10000), c=estimator.predict(data.reshape(-1,1)))
plt.show()




We see that the yellow dots are the 'normal' instances and the dark blue ones are the anomalies. You can also notice that this algorithm selects an anomaly threshold very similar to the one we obtained using classical statistics methods.

Non-parametric algorithms

This is a play on the Random Forests algorithm based on an ensemble of decision trees. The decision trees are trained to identify anomalies through information loss decisions and isolating instances of your data.

from sklearn.ensemble import IsolationForest
estimator = IsolationForest()
estimator.fit(data.reshape(-1,1))
plt.scatter(data, np.random.randint(0,1,10000), c=estimator.predict(data.reshape(-1,1)))
plt.show()




Non-parametric $p$-value estimation algorithms

These are personal favorite and have been very successful in a number of tasks.

Learning Minimum Volume Sets
http://www.stat.rice.edu/~cscott/pubs/minvol06jmlr.pdf

Anomaly Detection with Score functions based on Nearest Neighbor Graphs
https://arxiv.org/abs/0910.5461

New statistic in P-value estimation for anomaly detection
http://ieeexplore.ieee.org/document/6319713/



I don't have time tonight but tomorrow evening I will elaborate on the GLRT, and the $p$-value estimation algorithms.



Extending the above to more dimensions

Of course all the described algorithm can be extended to $n$-dimensions. This is hard to plot so we will go through a 2D example and a 3D example. However, determining the probability of these being anomalies is not so simple. For some of these algorithms you will need to devise some sort of distance metric.

2D example

Let's consider this artificial data. Where tall people above 170cm have eye colors with a blue content of mean 1 and variance 0.5. And short people have lower blue content at a mean 0 and variance 0.5 (let's just assume eyes can have negative pigmentation whatever that could possibly mean).

import numpy as np

n = 10000
X = np.zeros((n,2))
X[:,0] = np.random.normal(170, 15, n)
for i in range(n):
    if X[i,0] &lt; 170:
        X[i,1] = np.random.normal(0, 0.5)
    else:
        X[i,1] = np.random.normal(1, 0.5)

plt.scatter(X[:,0], X[:,1])
plt.show()




Classical statistics

Any p-value test used in classical statistics can be extended to multiple dimensions. This is the easiest when the features are uncorrelated. When this is the case the probability of each feature can be multiplied together. If this is not the case then you need to compute the p-value for a 2D Gaussian distribution. 

For our case this is even worse! It's bimodal (there are two means for a feature). You can see this here. Our data would look something like this image from the Wiki.



I hope you can appreciate how difficult this problem is to solve. Luckily this is no longer the classical era of statistics. 

Parametric techniques

These algorithms can easily be moved onto $n$-dimensions.

from sklearn.covariance import EllipticEnvelope
estimator = EllipticEnvelope()
estimator.fit(X)
plt.scatter(X[:,0], X[:,1], c=estimator.predict(X))
plt.show()




In the 1D example we discussed how this algorithm assumes a unimodal Gaussian distribution. This is not the case for our current data and you can see the failures of this model. It seems to classify many extreme examples as being part of the nominal set when they should have likely been identified as anomalies. 

Non-parametric techniques

This technique can also be applied to $n$-dimensional data. Moreover, it can be used for multimodal Gaussian distributions as you can see we get better results.

from sklearn.ensemble import IsolationForest
estimator = IsolationForest()
estimator.fit(X)
plt.scatter(X[:,0], X[:,1], c=estimator.predict(X))
plt.show()