use of de.lmu.ifi.dbs.elki.math.MeanVariance in project elki by elki-project.
the class GreedyEnsembleExperiment method run.
@Override
public void run() {
// Note: the database contains the *result vectors*, not the original data.
final Database database = inputstep.getDatabase();
Relation<NumberVector> relation = database.getRelation(TypeUtil.NUMBER_VECTOR_FIELD);
final Relation<String> labels = DatabaseUtil.guessLabelRepresentation(database);
final DBID firstid = DBIDUtil.deref(labels.iterDBIDs());
final String firstlabel = labels.get(firstid);
if (!firstlabel.matches("bylabel")) {
throw new AbortException("No 'by label' reference outlier found, which is needed for weighting!");
}
relation = applyPrescaling(prescaling, relation, firstid);
final int numcand = relation.size() - 1;
// Dimensionality and reference vector
final int dim = RelationUtil.dimensionality(relation);
final NumberVector refvec = relation.get(firstid);
// Build the positive index set for ROC AUC.
VectorNonZero positive = new VectorNonZero(refvec);
final int desired_outliers = (int) (rate * dim);
int union_outliers = 0;
final int[] outliers_seen = new int[dim];
// Merge the top-k for each ensemble member, until we have enough
// candidates.
{
int k = 0;
ArrayList<DecreasingVectorIter> iters = new ArrayList<>(numcand);
if (minvote >= numcand) {
minvote = Math.max(1, numcand - 1);
}
for (DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
// Skip "by label", obviously
if (DBIDUtil.equal(firstid, iditer)) {
continue;
}
iters.add(new DecreasingVectorIter(relation.get(iditer)));
}
loop: while (union_outliers < desired_outliers) {
for (DecreasingVectorIter iter : iters) {
if (!iter.valid()) {
LOG.warning("Union_outliers=" + union_outliers + " < desired_outliers=" + desired_outliers + " minvote=" + minvote);
break loop;
}
int cur = iter.dim();
outliers_seen[cur] += 1;
if (outliers_seen[cur] == minvote) {
union_outliers += 1;
}
iter.advance();
}
k++;
}
LOG.verbose("Merged top " + k + " outliers to: " + union_outliers + " outliers (desired: at least " + desired_outliers + ")");
}
// Build the final weight vector.
final double[] estimated_weights = new double[dim];
final double[] estimated_truth = new double[dim];
updateEstimations(outliers_seen, union_outliers, estimated_weights, estimated_truth);
DoubleVector estimated_truth_vec = DoubleVector.wrap(estimated_truth);
PrimitiveDistanceFunction<NumberVector> wdist = getDistanceFunction(estimated_weights);
PrimitiveDistanceFunction<NumberVector> tdist = wdist;
// Build the naive ensemble:
final double[] naiveensemble = new double[dim];
{
double[] buf = new double[numcand];
for (int d = 0; d < dim; d++) {
int i = 0;
for (DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
if (DBIDUtil.equal(firstid, iditer)) {
continue;
}
final NumberVector vec = relation.get(iditer);
buf[i] = vec.doubleValue(d);
i++;
}
naiveensemble[d] = voting.combine(buf, i);
if (Double.isNaN(naiveensemble[d])) {
LOG.warning("NaN after combining: " + FormatUtil.format(buf) + " i=" + i + " " + voting.toString());
}
}
}
DoubleVector naivevec = DoubleVector.wrap(naiveensemble);
// Compute single AUC scores and estimations.
// Remember the method most similar to the estimation
double bestauc = 0.0;
String bestaucstr = "";
double bestcost = Double.POSITIVE_INFINITY;
String bestcoststr = "";
DBID bestid = null;
double bestest = Double.POSITIVE_INFINITY;
{
final double[] greedyensemble = new double[dim];
// Compute individual scores
for (DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
if (DBIDUtil.equal(firstid, iditer)) {
continue;
}
// fout.append(labels.get(id));
final NumberVector vec = relation.get(iditer);
singleEnsemble(greedyensemble, vec);
double auc = ROCEvaluation.computeROCAUC(positive, new DecreasingVectorIter(DoubleVector.wrap(greedyensemble)));
double estimated = wdist.distance(DoubleVector.wrap(greedyensemble), estimated_truth_vec);
double cost = tdist.distance(DoubleVector.wrap(greedyensemble), refvec);
LOG.verbose("ROC AUC: " + auc + " estimated " + estimated + " cost " + cost + " " + labels.get(iditer));
if (auc > bestauc) {
bestauc = auc;
bestaucstr = labels.get(iditer);
}
if (cost < bestcost) {
bestcost = cost;
bestcoststr = labels.get(iditer);
}
if (estimated < bestest || bestid == null) {
bestest = estimated;
bestid = DBIDUtil.deref(iditer);
}
}
}
// Initialize ensemble with "best" method
if (prescaling != null) {
LOG.verbose("Input prescaling: " + prescaling);
}
LOG.verbose("Distance function: " + wdist);
LOG.verbose("Ensemble voting: " + voting);
if (scaling != null) {
LOG.verbose("Ensemble rescaling: " + scaling);
}
LOG.verbose("Initial estimation of outliers: " + union_outliers);
LOG.verbose("Initializing ensemble with: " + labels.get(bestid));
ModifiableDBIDs ensemble = DBIDUtil.newArray(bestid);
ModifiableDBIDs enscands = DBIDUtil.newHashSet(relation.getDBIDs());
ModifiableDBIDs dropped = DBIDUtil.newHashSet(relation.size());
dropped.add(firstid);
enscands.remove(bestid);
enscands.remove(firstid);
final double[] greedyensemble = new double[dim];
singleEnsemble(greedyensemble, relation.get(bestid));
// Greedily grow the ensemble
final double[] testensemble = new double[dim];
while (enscands.size() > 0) {
NumberVector greedyvec = DoubleVector.wrap(greedyensemble);
final double oldd = wdist.distance(estimated_truth_vec, greedyvec);
final int heapsize = enscands.size();
ModifiableDoubleDBIDList heap = DBIDUtil.newDistanceDBIDList(heapsize);
double[] tmp = new double[dim];
for (DBIDIter iter = enscands.iter(); iter.valid(); iter.advance()) {
final NumberVector vec = relation.get(iter);
singleEnsemble(tmp, vec);
double diversity = wdist.distance(DoubleVector.wrap(greedyensemble), greedyvec);
heap.add(diversity, iter);
}
heap.sort();
for (DoubleDBIDListMIter it = heap.iter(); heap.size() > 0; it.remove()) {
// Last
it.seek(heap.size() - 1);
enscands.remove(it);
final NumberVector vec = relation.get(it);
// Build combined ensemble.
{
double[] buf = new double[ensemble.size() + 1];
for (int i = 0; i < dim; i++) {
int j = 0;
for (DBIDIter iter = ensemble.iter(); iter.valid(); iter.advance()) {
buf[j] = relation.get(iter).doubleValue(i);
j++;
}
buf[j] = vec.doubleValue(i);
testensemble[i] = voting.combine(buf, j + 1);
}
}
applyScaling(testensemble, scaling);
NumberVector testvec = DoubleVector.wrap(testensemble);
double newd = wdist.distance(estimated_truth_vec, testvec);
// labels.get(bestadd));
if (newd < oldd) {
System.arraycopy(testensemble, 0, greedyensemble, 0, dim);
ensemble.add(it);
// Recompute heap
break;
} else {
dropped.add(it);
// logger.verbose("Discarding: " + labels.get(bestadd));
if (refine_truth) {
// Update target vectors and weights
ArrayList<DecreasingVectorIter> iters = new ArrayList<>(numcand);
for (DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
// Skip "by label", obviously
if (DBIDUtil.equal(firstid, iditer) || dropped.contains(iditer)) {
continue;
}
iters.add(new DecreasingVectorIter(relation.get(iditer)));
}
if (minvote >= iters.size()) {
minvote = iters.size() - 1;
}
union_outliers = 0;
Arrays.fill(outliers_seen, 0);
while (union_outliers < desired_outliers) {
for (DecreasingVectorIter iter : iters) {
if (!iter.valid()) {
break;
}
int cur = iter.dim();
if (outliers_seen[cur] == 0) {
outliers_seen[cur] = 1;
} else {
outliers_seen[cur] += 1;
}
if (outliers_seen[cur] == minvote) {
union_outliers += 1;
}
iter.advance();
}
}
LOG.warning("New num outliers: " + union_outliers);
updateEstimations(outliers_seen, union_outliers, estimated_weights, estimated_truth);
estimated_truth_vec = DoubleVector.wrap(estimated_truth);
}
}
}
}
// Build the improved ensemble:
StringBuilder greedylbl = new StringBuilder();
{
for (DBIDIter iter = ensemble.iter(); iter.valid(); iter.advance()) {
if (greedylbl.length() > 0) {
greedylbl.append(' ');
}
greedylbl.append(labels.get(iter));
}
}
DoubleVector greedyvec = DoubleVector.wrap(greedyensemble);
if (refine_truth) {
LOG.verbose("Estimated outliers remaining: " + union_outliers);
}
LOG.verbose("Greedy ensemble (" + ensemble.size() + "): " + greedylbl.toString());
LOG.verbose("Best single ROC AUC: " + bestauc + " (" + bestaucstr + ")");
LOG.verbose("Best single cost: " + bestcost + " (" + bestcoststr + ")");
// Evaluate the naive ensemble and the "shrunk" ensemble
double naiveauc, naivecost;
{
naiveauc = ROCEvaluation.computeROCAUC(positive, new DecreasingVectorIter(naivevec));
naivecost = tdist.distance(naivevec, refvec);
LOG.verbose("Naive ensemble AUC: " + naiveauc + " cost: " + naivecost);
LOG.verbose("Naive ensemble Gain: " + gain(naiveauc, bestauc, 1) + " cost gain: " + gain(naivecost, bestcost, 0));
}
double greedyauc, greedycost;
{
greedyauc = ROCEvaluation.computeROCAUC(positive, new DecreasingVectorIter(greedyvec));
greedycost = tdist.distance(greedyvec, refvec);
LOG.verbose("Greedy ensemble AUC: " + greedyauc + " cost: " + greedycost);
LOG.verbose("Greedy ensemble Gain to best: " + gain(greedyauc, bestauc, 1) + " cost gain: " + gain(greedycost, bestcost, 0));
LOG.verbose("Greedy ensemble Gain to naive: " + gain(greedyauc, naiveauc, 1) + " cost gain: " + gain(greedycost, naivecost, 0));
}
{
MeanVariance meanauc = new MeanVariance();
MeanVariance meancost = new MeanVariance();
HashSetModifiableDBIDs candidates = DBIDUtil.newHashSet(relation.getDBIDs());
candidates.remove(firstid);
for (int i = 0; i < 1000; i++) {
// Build the improved ensemble:
final double[] randomensemble = new double[dim];
{
DBIDs random = DBIDUtil.randomSample(candidates, ensemble.size(), (long) i);
double[] buf = new double[random.size()];
for (int d = 0; d < dim; d++) {
int j = 0;
for (DBIDIter iter = random.iter(); iter.valid(); iter.advance()) {
assert (!DBIDUtil.equal(firstid, iter));
final NumberVector vec = relation.get(iter);
buf[j] = vec.doubleValue(d);
j++;
}
randomensemble[d] = voting.combine(buf, j);
}
}
applyScaling(randomensemble, scaling);
NumberVector randomvec = DoubleVector.wrap(randomensemble);
double auc = ROCEvaluation.computeROCAUC(positive, new DecreasingVectorIter(randomvec));
meanauc.put(auc);
double cost = tdist.distance(randomvec, refvec);
meancost.put(cost);
}
LOG.verbose("Random ensemble AUC: " + meanauc.getMean() + " + stddev: " + meanauc.getSampleStddev() + " = " + (meanauc.getMean() + meanauc.getSampleStddev()));
LOG.verbose("Random ensemble Gain: " + gain(meanauc.getMean(), bestauc, 1));
LOG.verbose("Greedy improvement: " + (greedyauc - meanauc.getMean()) / meanauc.getSampleStddev() + " standard deviations.");
LOG.verbose("Random ensemble Cost: " + meancost.getMean() + " + stddev: " + meancost.getSampleStddev() + " = " + (meancost.getMean() + meanauc.getSampleStddev()));
LOG.verbose("Random ensemble Gain: " + gain(meancost.getMean(), bestcost, 0));
LOG.verbose("Greedy improvement: " + (meancost.getMean() - greedycost) / meancost.getSampleStddev() + " standard deviations.");
LOG.verbose("Naive ensemble Gain to random: " + gain(naiveauc, meanauc.getMean(), 1) + " cost gain: " + gain(naivecost, meancost.getMean(), 0));
LOG.verbose("Random ensemble Gain to naive: " + gain(meanauc.getMean(), naiveauc, 1) + " cost gain: " + gain(meancost.getMean(), naivecost, 0));
LOG.verbose("Greedy ensemble Gain to random: " + gain(greedyauc, meanauc.getMean(), 1) + " cost gain: " + gain(greedycost, meancost.getMean(), 0));
}
}
use of de.lmu.ifi.dbs.elki.math.MeanVariance in project elki by elki-project.
the class IntrinsicNearestNeighborAffinityMatrixBuilder method computePij.
/**
* Compute the sparse pij using the nearest neighbors only.
*
* @param ids ID range
* @param knnq kNN query
* @param square Use squared distances
* @param numberOfNeighbours Number of neighbors to get
* @param pij Output of distances
* @param indices Output of indexes
* @param initialScale Initial scaling factor
*/
protected void computePij(DBIDRange ids, KNNQuery<?> knnq, boolean square, int numberOfNeighbours, double[][] pij, int[][] indices, double initialScale) {
Duration timer = LOG.isStatistics() ? LOG.newDuration(this.getClass().getName() + ".runtime.neighborspijmatrix").begin() : null;
final double logPerp = FastMath.log(perplexity);
// Scratch arrays, resizable
DoubleArray dists = new DoubleArray(numberOfNeighbours + 10);
IntegerArray inds = new IntegerArray(numberOfNeighbours + 10);
// Compute nearest-neighbor sparse affinity matrix
FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Finding neighbors and optimizing perplexity", ids.size(), LOG) : null;
MeanVariance mv = LOG.isStatistics() ? new MeanVariance() : null;
Mean mid = LOG.isStatistics() ? new Mean() : null;
for (DBIDArrayIter ix = ids.iter(); ix.valid(); ix.advance()) {
dists.clear();
inds.clear();
KNNList neighbours = knnq.getKNNForDBID(ix, numberOfNeighbours + 1);
convertNeighbors(ids, ix, square, neighbours, dists, inds, mid);
double beta = computeSigma(//
ix.getOffset(), //
dists, //
perplexity, //
logPerp, pij[ix.getOffset()] = new double[dists.size()]);
if (mv != null) {
// Sigma
mv.put(beta > 0 ? FastMath.sqrt(.5 / beta) : 0.);
}
indices[ix.getOffset()] = inds.toArray();
LOG.incrementProcessed(prog);
}
LOG.ensureCompleted(prog);
if (mid != null) {
LOG.statistics(new DoubleStatistic(getClass() + ".average-original-id", mid.getMean()));
}
// Sum of the sparse affinity matrix:
double sum = 0.;
for (int i = 0; i < pij.length; i++) {
final double[] pij_i = pij[i];
for (int offi = 0; offi < pij_i.length; offi++) {
int j = indices[i][offi];
if (j > i) {
// Exploit symmetry.
continue;
}
assert (i != j);
int offj = containsIndex(indices[j], i);
if (offj >= 0) {
// Found
sum += FastMath.sqrt(pij_i[offi] * pij[j][offj]);
}
}
}
final double scale = initialScale / (2 * sum);
for (int i = 0; i < pij.length; i++) {
final double[] pij_i = pij[i];
for (int offi = 0; offi < pij_i.length; offi++) {
int j = indices[i][offi];
assert (i != j);
int offj = containsIndex(indices[j], i);
if (offj >= 0) {
// Found
assert (indices[j][offj] == i);
// Exploit symmetry:
if (i < j) {
// Symmetrize
final double val = FastMath.sqrt(pij_i[offi] * pij[j][offj]);
pij_i[offi] = pij[j][offj] = MathUtil.max(val * scale, MIN_PIJ);
}
} else {
// Not found, so zero.
pij_i[offi] = 0;
}
}
}
if (LOG.isStatistics()) {
// timer != null, mv != null
LOG.statistics(timer.end());
LOG.statistics(new DoubleStatistic(NearestNeighborAffinityMatrixBuilder.class.getName() + ".sigma.average", mv.getMean()));
LOG.statistics(new DoubleStatistic(NearestNeighborAffinityMatrixBuilder.class.getName() + ".sigma.stddev", mv.getSampleStddev()));
}
}
use of de.lmu.ifi.dbs.elki.math.MeanVariance in project elki by elki-project.
the class EvaluateRankingQuality method run.
@Override
public HistogramResult run(Database database) {
final Relation<V> relation = database.getRelation(getInputTypeRestriction()[0]);
final DistanceQuery<V> distQuery = database.getDistanceQuery(relation, getDistanceFunction());
final KNNQuery<V> knnQuery = database.getKNNQuery(distQuery, relation.size());
if (LOG.isVerbose()) {
LOG.verbose("Preprocessing clusters...");
}
// Cluster by labels
Collection<Cluster<Model>> split = (new ByLabelOrAllInOneClustering()).run(database).getAllClusters();
// Compute cluster averages and covariance matrix
HashMap<Cluster<?>, double[]> averages = new HashMap<>(split.size());
HashMap<Cluster<?>, double[][]> covmats = new HashMap<>(split.size());
for (Cluster<?> clus : split) {
CovarianceMatrix covmat = CovarianceMatrix.make(relation, clus.getIDs());
averages.put(clus, covmat.getMeanVector());
covmats.put(clus, covmat.destroyToPopulationMatrix());
}
MeanVarianceStaticHistogram hist = new MeanVarianceStaticHistogram(numbins, 0.0, 1.0);
if (LOG.isVerbose()) {
LOG.verbose("Processing points...");
}
FiniteProgress rocloop = LOG.isVerbose() ? new FiniteProgress("Computing ROC AUC values", relation.size(), LOG) : null;
ROCEvaluation roc = new ROCEvaluation();
// sort neighbors
for (Cluster<?> clus : split) {
ModifiableDoubleDBIDList cmem = DBIDUtil.newDistanceDBIDList(clus.size());
double[] av = averages.get(clus);
double[][] covm = covmats.get(clus);
for (DBIDIter iter = clus.getIDs().iter(); iter.valid(); iter.advance()) {
double d = mahalanobisDistance(covm, relation.get(iter).toArray(), av);
cmem.add(d, iter);
}
cmem.sort();
for (DBIDArrayIter it = cmem.iter(); it.valid(); it.advance()) {
KNNList knn = knnQuery.getKNNForDBID(it, relation.size());
double result = EvaluateClustering.evaluateRanking(roc, clus, knn);
hist.put(((double) it.getOffset()) / clus.size(), result);
LOG.incrementProcessed(rocloop);
}
}
LOG.ensureCompleted(rocloop);
// Collections.sort(results);
// Transform Histogram into a Double Vector array.
Collection<double[]> res = new ArrayList<>(relation.size());
for (ObjHistogram.Iter<MeanVariance> iter = hist.iter(); iter.valid(); iter.advance()) {
res.add(new double[] { iter.getCenter(), iter.getValue().getCount(), iter.getValue().getMean(), iter.getValue().getSampleVariance() });
}
return new HistogramResult("Ranking Quality Histogram", "ranking-histogram", res);
}
use of de.lmu.ifi.dbs.elki.math.MeanVariance in project elki by elki-project.
the class RangeQuerySelectivity method run.
public Result run(Database database, Relation<V> relation) {
DistanceQuery<V> distQuery = database.getDistanceQuery(relation, getDistanceFunction());
RangeQuery<V> rangeQuery = database.getRangeQuery(distQuery, radius);
MeanVariance numres = new MeanVariance();
final DBIDs ids = DBIDUtil.randomSample(relation.getDBIDs(), sampling, random);
FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Performing range queries", ids.size(), LOG) : null;
for (DBIDIter iter = ids.iter(); iter.valid(); iter.advance()) {
numres.put(rangeQuery.getRangeForDBID(iter, radius).size());
LOG.incrementProcessed(prog);
}
LOG.ensureCompleted(prog);
final String prefix = this.getClass().getName();
LOG.statistics(new DoubleStatistic(prefix + ".mean", numres.getMean()));
LOG.statistics(new DoubleStatistic(prefix + ".std", numres.getSampleStddev()));
LOG.statistics(new DoubleStatistic(prefix + ".norm.mean", numres.getMean() / relation.size()));
LOG.statistics(new DoubleStatistic(prefix + ".norm.std", numres.getSampleStddev() / relation.size()));
LOG.statistics(new LongStatistic(prefix + ".samplesize", ids.size()));
return null;
}
use of de.lmu.ifi.dbs.elki.math.MeanVariance in project elki by elki-project.
the class KDEOS method computeOutlierScores.
/**
* Compute the final KDEOS scores.
*
* @param knnq kNN query
* @param ids IDs to process
* @param densities Density estimates
* @param kdeos Score outputs
* @param minmax Minimum and maximum scores
*/
protected void computeOutlierScores(KNNQuery<O> knnq, final DBIDs ids, WritableDataStore<double[]> densities, WritableDoubleDataStore kdeos, DoubleMinMax minmax) {
final int knum = kmax + 1 - kmin;
FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Computing KDEOS scores", ids.size(), LOG) : null;
double[][] scratch = new double[knum][kmax + 5];
MeanVariance mv = new MeanVariance();
for (DBIDIter iter = ids.iter(); iter.valid(); iter.advance()) {
double[] dens = densities.get(iter);
KNNList neighbors = knnq.getKNNForDBID(iter, kmax + 1);
if (scratch[0].length < neighbors.size()) {
// Resize scratch. Add some extra margin again.
scratch = new double[knum][neighbors.size() + 5];
}
{
// Store density matrix of neighbors
int i = 0;
for (DoubleDBIDListIter neighbor = neighbors.iter(); neighbor.valid(); neighbor.advance(), i++) {
double[] ndens = densities.get(neighbor);
for (int k = 0; k < knum; k++) {
scratch[k][i] = ndens[k];
}
}
assert (i == neighbors.size());
}
// Compute means and stddevs for each k
double score = 0.;
for (int i = 0; i < knum; i++) {
mv.reset();
for (int j = 0; j < neighbors.size(); j++) {
mv.put(scratch[i][j]);
}
final double mean = mv.getMean(), stddev = mv.getSampleStddev();
if (stddev > 0.) {
score += (mean - dens[i]) / stddev;
}
}
// average
score /= knum;
score = NormalDistribution.standardNormalCDF(score);
minmax.put(score);
kdeos.put(iter, score);
LOG.incrementProcessed(prog);
}
LOG.ensureCompleted(prog);
}
Aggregations