From 9325f04856b8cd90dfc504010ae08f4c9c194bf1 Mon Sep 17 00:00:00 2001
From: Daniel Waggoner <dwaggoner@frbatlanta.org>
Date: Thu, 3 May 2012 17:29:24 -0400
Subject: [PATCH] changes to elliptical testing and minor changes of
dw_elliptical.c (cherry picked from commit
f977618b9da25a83daaf1bc1eb09cf6a71cd2c32)
---
elliptical/dw_elliptical.c | 276 +++-----------------------------
elliptical/dw_elliptical_test.c | 236 +++++++++++++++++++++++++++
include/dw_elliptical.h | 7 +-
3 files changed, 260 insertions(+), 259 deletions(-)
create mode 100644 elliptical/dw_elliptical_test.c
diff --git a/elliptical/dw_elliptical.c b/elliptical/dw_elliptical.c
index 11cfd66..89a94f6 100644
--- a/elliptical/dw_elliptical.c
+++ b/elliptical/dw_elliptical.c
@@ -170,7 +170,7 @@ static void print_info_truncated_gaussian(FILE *f_out, TElliptical *elliptical)
}
/*
- Creates a power distribution on the elliptical annulus which is given by the
+ Creates a truncated gaussian distribution on the elliptical annulus which is given by the
set of all z such that
r = sqrt((z - center)' * Inverse(scale * scale') * (z - center))
@@ -522,6 +522,16 @@ static PRECISION logdensity_radius_step(PRECISION radius, TElliptical *elliptica
}
}
+static PRECISION cummulative_radius_step(PRECISION radius, TElliptical *elliptical)
+{
+ TElliptical_step *d=(TElliptical_step*)elliptical->data;
+ int i;
+ if (radius <= d->t[0]) return 0.0;
+ for (i=1; i <= d->m; i++)
+ if (radius <= d->t[i]) return ((PRECISION)(i-1) + (radius - d->t[i-1])/(d->t[i] - d->t[i-1]))/(PRECISION)d->m;
+ return 1.0;
+}
+
static PRECISION draw_step(PRECISION *draw, TElliptical *elliptical)
{
PRECISION r, s;
@@ -543,8 +553,15 @@ static PRECISION draw_step(PRECISION *draw, TElliptical *elliptical)
static void print_info_step(FILE *f_out, TElliptical *elliptical)
{
TElliptical_step *d=(TElliptical_step*)(elliptical->data);
- if (f_out) fprintf(f_out,"%s\n dim=%d\n lower bound=%lg\n upper bound=%lg\n",
- elliptical->type,elliptical->dim,d->t[0],d->t[d->m]);
+ int i;
+ if (f_out)
+ {
+ fprintf(f_out,"%s\n dim=%d\n lower bound=%lg\n upper bound=%lg\n",elliptical->type,elliptical->dim,d->t[0],d->t[d->m]);
+ fprintf(f_out,"table size=%d\n table=\n",d->m);
+ for (i=0; i <= d->m; i++) fprintf(f_out,"%lg ",d->t[i]);
+ fprintf(f_out,"\n");
+ }
+
}
/*
@@ -597,6 +614,7 @@ TElliptical* CreateElliptical_Step(int dim, TVector center, TMatrix scale, PRECI
elliptical->data=(void*)data;
elliptical->logdensity_radius=logdensity_radius_step;
+ elliptical->cummulative_radius=cummulative_radius_step;
elliptical->logdensity_draw=logdensity_draw;
elliptical->draw_vector=draw_step;
elliptical->print_info=print_info_step;
@@ -608,256 +626,4 @@ TElliptical* CreateElliptical_Step(int dim, TVector center, TMatrix scale, PRECI
/*******************************************************************************/
/*******************************************************************************/
/*******************************************************************************/
-
-
-/*******************************************************************************/
-/******************************** Test Routines ********************************/
-/*******************************************************************************
-#include "dw_matrix_sort.h"
-#include <stdio.h>
-static void TestSimulateEllipticalCummulativeDensity(TElliptical *elliptical, int count, char *filename)
-{
- TVector R;
- int d=200, i, j;
- FILE *f_out;
- PRECISION min, max, inc, x;
-
- if (!filename) filename="tmp.csv";
- if (elliptical->cummulative_radius && (f_out=fopen(filename,"wt")))
- {
- R=CreateVector(count);
- for (i=0; i < count; i++)
- ElementV(R,i)=DrawElliptical((PRECISION*)NULL,elliptical);
- SortVectorAscending(R,R);
-
- min=ElementV(R,0);
- max=ElementV(R,count-1);
- inc=(max-min)/(PRECISION)(d-1);
-
- for (x=min, i=j=0; i < d; x+=inc, i++)
- {
- for ( ; j < count; j++)
- if (ElementV(R,j) >= x) break;
- fprintf(f_out,"%lf,%lf,%lf\n",x,(double)j/(double)count,CummulativeDensityElliptical_Radius(x,elliptical));
- }
- fclose(f_out);
- }
-}
-
-// Prints, in a format suitable for graphing, the cummulative density.
-void TestEllipticalCummulativeDensity(void)
-{
- TVector center;
- TMatrix scale, variance;
- PRECISION lower_bound=2.0, upper_bound=10.0, power=5, s=2.3, *table;
- int i, j, dim=3, count=50000, m=10;
- char *filename="cummulative.csv";
- TElliptical *elliptical;
-
- dw_initialize_generator(0);
-
- dw_NormalVector(center=CreateVector(dim));
- dw_NormalMatrix(scale=CreateMatrix(dim,dim));
-
- printf("Cummulative density test\n");
- printf("Count = %d\n",count);
-
- // Truncated gaussian
- variance=ProductTransposeMM((TMatrix)NULL,scale,scale);
- elliptical=CreateElliptical_TruncatedGaussian(dim,center,variance,lower_bound,upper_bound,s);
- printf("\nCumulative density test of truncated gaussian distribution...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalCummulativeDensity(elliptical,count,filename);
- FreeMatrix(variance);
- FreeElliptical(elliptical);
-
- // Gaussian
- variance=ProductTransposeMM((TMatrix)NULL,scale,scale);
- elliptical=CreateElliptical_Gaussian(dim,center,variance);
- printf("\nCumulative density test of gaussian distribution...\n");
- PrintEllipticalInfo(stdout,elliptical);
- printf("\nGaussian test...\n");
- TestSimulateEllipticalCummulativeDensity(elliptical,count,filename);
- FreeMatrix(variance);
- FreeElliptical(elliptical);
-
- // Truncated Power
- elliptical=CreateElliptical_TruncatedPower(dim,center,scale,lower_bound,upper_bound,power);
- printf("\nCumulative density test of truncated power distribution...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalCummulativeDensity(elliptical,count,filename);
- FreeElliptical(elliptical);
-
- // Uniform
- elliptical=CreateElliptical_Uniform(dim,center,scale,lower_bound,upper_bound);
- printf("\nCumulative density test of uniform distribution...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalCummulativeDensity(elliptical,count,filename);
- FreeElliptical(elliptical);
-
- // Power
- elliptical=CreateElliptical_Power(dim,center,scale,upper_bound,power);
- printf("\nCumulative density test of power distribution...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalCummulativeDensity(elliptical,count,filename);
- FreeElliptical(elliptical);
-
- // Step
- table=(PRECISION*)dw_malloc((m+1)*sizeof(PRECISION));
- table[0]=lower_bound;
- for (i=1; i < m; i++)
- table[i]=table[i-1]+dw_uniform_rnd()*(upper_bound-table[i-1]);
- table[m]=upper_bound;
- factor=sqrt(VarianceScale_TruncatedPower(dim,lower_bound,upper_bound,power));
- ProductMS(scale,scale,1.0/factor);
- elliptical=CreateElliptical_Step(dim,center,scale,table,m);
- ProductMS(scale,scale,factor);
- printf("\nCumulative density test of step distribution...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- dw_free(table);
- FreeElliptical(elliptical);
-
- FreeMatrix(scale);
- FreeVector(center);
-}
-
-static void TestSimulateEllipticalVariance(TElliptical *elliptical, int count)
-{
- TVector mean, v;
- TMatrix variance, x;
- int i;
-
- // Create matrices and vectors
- InitializeVector(mean=CreateVector(elliptical->dim),0.0);
- InitializeMatrix(variance=CreateMatrix(elliptical->dim,elliptical->dim),0.0);
- v=CreateVector(elliptical->dim);
- x=CreateMatrix(elliptical->dim,elliptical->dim);
-
- // Simulate mean and variance
- for (i=0; i < count; i++)
- {
- DrawElliptical(pElementV(v),elliptical);
- AddVV(mean,mean,v);
- OuterProduct(x,v,v);
- AddMM(variance,variance,x);
- }
- ProductVS(mean,mean,1/(PRECISION)count);
- ProductMS(variance,variance,1/(PRECISION)count);
- OuterProduct(x,mean,mean);
- SubtractMM(variance,variance,x);
-
- // Difference in centers
- SubtractVV(v,elliptical->center,mean);
- printf("difference in estimated vs acutal center: %lf\n",Norm(v));
-
- // Difference in variances
- ProductMS(x,elliptical->scale,sqrt(elliptical->variance_scale));
- ProductTransposeMM(x,x,x);
- SubtractMM(x,x,variance);
- printf("difference in estimated vs actual variance: %lf\n",MatrixNorm(x));
-
- // Prints the variance matrices
- ProductMS(x,elliptical->scale,sqrt(elliptical->variance_scale));
- ProductTransposeMM(x,x,x);
- printf("Actual Variance:\n");
- dw_PrintMatrix(stdout,x,"%lg ");
- printf("Estimated Variance:\n");
- dw_PrintMatrix(stdout,variance,"%lg ");
-
- // Clean up
- FreeMatrix(variance);
- FreeVector(mean);
- FreeVector(v);
- FreeMatrix(x);
-}
-
-// Prints the difference difference, in terms of norms, between estimated and
-// actual variances and means.
-void TestEllipticalVariance(void)
-{
- TMatrix scale, variance;
- TVector center;
- PRECISION factor, lower_bound=0.0, upper_bound=10.0, power=5.0, s=2.0, *table;
- int i, j, dim=2, count=50000, m=10;
- TElliptical *elliptical;
-
- dw_initialize_generator(0);
-
- dw_NormalVector(center=CreateVector(dim));
- dw_NormalMatrix(scale=CreateMatrix(dim,dim));
- variance=ProductTransposeMM((TMatrix)NULL,scale,scale);
-
- printf("Mean and Variance test\n");
- printf("Count = %d\n",count);
-
- // Gaussian family
- factor=VarianceScale_TruncatedGaussian(dim,lower_bound,upper_bound,s);
- ProductMS(variance,variance,1.0/factor);
- elliptical=CreateElliptical_TruncatedGaussian(dim,center,variance,lower_bound,upper_bound,s);
- ProductMS(variance,variance,factor);
- printf("\nTruncated gaussian test (factor=%lf)...\n",factor);
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- FreeElliptical(elliptical);
-
- factor=sqrt(VarianceScale_TruncatedGaussian(dim,0.0,PLUS_INFINITY,1.0));
- elliptical=CreateElliptical_Gaussian(dim,center,variance);
- printf("\nGaussian test...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- FreeElliptical(elliptical);
-
- // Truncated Power family
- factor=sqrt(VarianceScale_TruncatedPower(dim,lower_bound,upper_bound,power));
- ProductMS(scale,scale,1.0/factor);
- elliptical=CreateElliptical_TruncatedPower(dim,center,scale,lower_bound,upper_bound,power);
- ProductMS(scale,scale,factor);
- printf("\nTruncated power test...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- FreeElliptical(elliptical);
-
- factor=sqrt(VarianceScale_TruncatedPower(dim,lower_bound,upper_bound,dim));
- ProductMS(scale,scale,1.0/factor);
- elliptical=CreateElliptical_Uniform(dim,center,scale,lower_bound,upper_bound);
- ProductMS(scale,scale,factor);
- printf("\nUniform test...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- FreeElliptical(elliptical);
-
- factor=sqrt(VarianceScale_TruncatedPower(dim,0.0,upper_bound,power));
- ProductMS(scale,scale,1.0/factor);
- elliptical=CreateElliptical_Power(dim,center,scale,upper_bound,power);
- ProductMS(scale,scale,factor);
- printf("\nPower test...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- FreeElliptical(elliptical);
-
- // Step
- table=(PRECISION*)dw_malloc((m+1)*sizeof(PRECISION));
- table[0]=lower_bound;
- for (i=1; i < m; i++)
- table[i]=table[i-1]+dw_uniform_rnd()*(upper_bound-table[i-1]);
- table[m]=upper_bound;
- factor=sqrt(VarianceScale_Step(dim,table,m));
- ProductMS(scale,scale,1.0/factor);
- elliptical=CreateElliptical_Step(dim,center,scale,table,m);
- ProductMS(scale,scale,factor);
- printf("\nStep test...\n");
- PrintEllipticalInfo(stdout,elliptical);
- TestSimulateEllipticalVariance(elliptical,count);
- dw_free(table);
- FreeElliptical(elliptical);
-
- // Clean up
- FreeMatrix(variance);
- FreeVector(center);
- FreeMatrix(scale);
-}
-*******************************************************************************/
-/*******************************************************************************/
-/*******************************************************************************/
#undef PI
diff --git a/elliptical/dw_elliptical_test.c b/elliptical/dw_elliptical_test.c
new file mode 100644
index 0000000..9f93319
--- /dev/null
+++ b/elliptical/dw_elliptical_test.c
@@ -0,0 +1,236 @@
+/*******************************************************************************/
+/************************** Elliptical Test Routines ***************************/
+/*******************************************************************************/
+#include "dw_matrix_sort.h"
+#include "dw_elliptical.h"
+#include "dw_rand.h"
+#include "dw_matrix_rand.h"
+#include "dw_parse_cmd.h"
+#include "dw_std.h"
+
+#include <stdio.h>
+#include <math.h>
+
+static void TestSimulateEllipticalVariance(TElliptical *elliptical, int count, FILE *f_out);
+static void TestSimulateEllipticalCummulativeDensity(TElliptical *elliptical, int count, FILE *f_out);
+
+static void TestSimulateEllipticalCummulativeDensity(TElliptical *elliptical, int count, FILE *f_out)
+{
+ TVector R;
+ int d=200, i, j;
+ PRECISION min, max, inc, x;
+
+ if (!f_out)
+ {
+ printf("Null file pointer\n");
+ exit(0);
+ }
+
+ PrintEllipticalInfo(f_out,elliptical);
+
+ if (elliptical->cummulative_radius)
+ {
+ R=CreateVector(count);
+ for (i=0; i < count; i++)
+ ElementV(R,i)=DrawElliptical((PRECISION*)NULL,elliptical);
+ SortVectorAscending(R,R);
+
+ min=ElementV(R,0);
+ max=ElementV(R,count-1);
+ inc=(max-min)/(PRECISION)(d-1);
+
+ for (x=min, i=j=0; i < d; x+=inc, i++)
+ {
+ for ( ; j < count; j++)
+ if (ElementV(R,j) >= x) break;
+ fprintf(f_out,"%lf,%lf,%lf\n",x,(double)j/(double)count,CummulativeDensityElliptical_Radius(x,elliptical));
+ }
+ }
+
+ fprintf(f_out,"\n\n");
+}
+
+static void TestSimulateEllipticalVariance(TElliptical *elliptical, int count, FILE *f_out)
+{
+ TVector mean, v;
+ TMatrix variance, x;
+ int i;
+
+ if (!f_out)
+ {
+ printf("Null file pointer\n");
+ exit(0);
+ }
+
+ PrintEllipticalInfo(f_out,elliptical);
+
+ // Create matrices and vectors
+ InitializeVector(mean=CreateVector(elliptical->dim),0.0);
+ InitializeMatrix(variance=CreateMatrix(elliptical->dim,elliptical->dim),0.0);
+ v=CreateVector(elliptical->dim);
+ x=CreateMatrix(elliptical->dim,elliptical->dim);
+
+ // Simulate mean and variance
+ for (i=0; i < count; i++)
+ {
+ DrawElliptical(pElementV(v),elliptical);
+ AddVV(mean,mean,v);
+ OuterProduct(x,v,v);
+ AddMM(variance,variance,x);
+ }
+ ProductVS(mean,mean,1/(PRECISION)count);
+ ProductMS(variance,variance,1/(PRECISION)count);
+ OuterProduct(x,mean,mean);
+ SubtractMM(variance,variance,x);
+
+ // Difference in centers
+ fprintf(f_out,"acutal mean\n");
+ dw_PrintVector(f_out,elliptical->center,"%lf ");
+ fprintf(f_out,"difference between estimated and actual mean divided by standard deviation\n");
+ SubtractVV(v,elliptical->center,mean);
+ for (i=0; i < elliptical->dim; i++)
+ fprintf(f_out,"%8.5lf ",ElementV(v,i)/sqrt(ElementM(variance,i,i)/count));
+ fprintf(f_out,"\n");
+
+ // Difference in variances
+ ProductMS(x,elliptical->scale,sqrt(elliptical->variance_scale));
+ ProductTransposeMM(x,x,x);
+ SubtractMM(x,x,variance);
+ fprintf(f_out,"norm of the difference between estimated and actual variance: %lf\n",MatrixNorm(x));
+
+ // Prints the variance matrices
+ fprintf(f_out,"Actual Variance:\n");
+ ProductMS(x,elliptical->scale,sqrt(elliptical->variance_scale));
+ ProductTransposeMM(x,x,x);
+ dw_PrintMatrix(f_out,x,"%lg ");
+ fprintf(f_out,"difference between estimated and actual variance:\n");
+ SubtractMM(x,x,variance);
+ dw_PrintMatrix(f_out,x,"%lg ");
+
+ fprintf(f_out,"\n\n");
+
+ // Clean up
+ FreeMatrix(variance);
+ FreeVector(mean);
+ FreeVector(v);
+ FreeMatrix(x);
+}
+
+// Prints the difference, in terms of norms, between estimated and
+// actual variances and means.
+int main(int nargs, char **args)
+{
+ TMatrix base_scale, scale, variance;
+ TVector center;
+ PRECISION lower_bound, upper_bound, power, s, *table;
+ int i, seed, dim, count, table_size, cummulative;
+ TElliptical *elliptical;
+ FILE *f_out;
+
+ count=dw_ParseInteger_String(nargs,args,"ndraws",10000);
+ dim=dw_ParseInteger_String(nargs,args,"dim",3);
+ lower_bound=dw_ParseFloating_String(nargs,args,"lower_bound",2.0);
+ if (lower_bound < 0.0) lower_bound=0.0;
+ upper_bound=dw_ParseFloating_String(nargs,args,"upper_bound",8.0);
+ if (lower_bound >= upper_bound) upper_bound=lower_bound+1.0;
+ power=dw_ParseFloating_String(nargs,args,"power",5.0);
+ table_size=dw_ParseInteger_String(nargs,args,"table_size",10);
+ s=dw_ParseFloating_String(nargs,args,"s",1.0);
+ cummulative=(dw_FindArgument_String(nargs,args,"cummulative") == -1) ? 0 : 1;
+ seed=dw_ParseInteger_String(nargs,args,"seed",0);
+
+ f_out=cummulative ? fopen("cummulative.csv","wt") : stdout;
+
+ dw_initialize_generator(seed);
+
+ dw_NormalVector(center=CreateVector(dim));
+ dw_NormalMatrix(base_scale=CreateMatrix(dim,dim));
+ variance=CreateMatrix(dim,dim);
+ scale=CreateMatrix(dim,dim);
+
+ fprintf(f_out,"Mean and Variance test\n");
+ fprintf(f_out,"Count = %d\n",count);
+
+ // Gaussian family
+ variance=ProductTransposeMM((TMatrix)NULL,base_scale,base_scale);
+ ProductMS(variance,variance,1.0/VarianceScale_TruncatedGaussian(dim,lower_bound,upper_bound,s));
+ elliptical=CreateElliptical_TruncatedGaussian(dim,center,variance,lower_bound,upper_bound,s);
+ if (cummulative)
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ else
+ TestSimulateEllipticalVariance(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ //==========
+ if (cummulative)
+ {
+ variance=ProductTransposeMM((TMatrix)NULL,base_scale,base_scale);
+ elliptical=CreateElliptical_TruncatedGaussian(dim,center,variance,lower_bound,upper_bound,s);
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ elliptical=CreateElliptical_TruncatedGaussian(dim,center,variance,lower_bound,upper_bound,s+100);
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+ }
+ //============
+
+ variance=ProductTransposeMM((TMatrix)NULL,base_scale,base_scale);
+ elliptical=CreateElliptical_Gaussian(dim,center,variance);
+ if (cummulative)
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ else
+ TestSimulateEllipticalVariance(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ // Truncated Power family
+ ProductMS(scale,base_scale,1.0/sqrt(VarianceScale_TruncatedPower(dim,lower_bound,upper_bound,power)));
+ elliptical=CreateElliptical_TruncatedPower(dim,center,scale,lower_bound,upper_bound,power);
+ if (cummulative)
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ else
+ TestSimulateEllipticalVariance(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ ProductMS(scale,base_scale,1.0/sqrt(VarianceScale_TruncatedPower(dim,lower_bound,upper_bound,dim)));
+ elliptical=CreateElliptical_Uniform(dim,center,scale,lower_bound,upper_bound);
+ if (cummulative)
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ else
+ TestSimulateEllipticalVariance(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ ProductMS(scale,base_scale,1.0/sqrt(VarianceScale_TruncatedPower(dim,0.0,upper_bound,power)));
+ elliptical=CreateElliptical_Power(dim,center,scale,upper_bound,power);
+ if (cummulative)
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ else
+ TestSimulateEllipticalVariance(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ // Step
+ table=(PRECISION*)dw_malloc((table_size+1)*sizeof(PRECISION));
+ table[0]=lower_bound;
+ for (i=1; i < table_size; i++)
+ table[i]=table[i-1]+dw_uniform_rnd()*(upper_bound-table[i-1]);
+ table[table_size]=upper_bound;
+ ProductMS(scale,base_scale,1.0/sqrt(VarianceScale_Step(dim,table,table_size)));
+ elliptical=CreateElliptical_Step(dim,center,scale,table,table_size);
+ dw_free(table);
+ if (cummulative)
+ TestSimulateEllipticalCummulativeDensity(elliptical,count,f_out);
+ else
+ TestSimulateEllipticalVariance(elliptical,count,f_out);
+ FreeElliptical(elliptical);
+
+ // Clean up
+ FreeMatrix(variance);
+ FreeVector(center);
+ FreeMatrix(scale);
+ FreeMatrix(base_scale);
+
+ return 0;
+}
+/*******************************************************************************/
+/*******************************************************************************/
+/*******************************************************************************/
diff --git a/include/dw_elliptical.h b/include/dw_elliptical.h
index d6c3b7e..8d0f595 100644
--- a/include/dw_elliptical.h
+++ b/include/dw_elliptical.h
@@ -84,6 +84,7 @@ void PrintEllipticalInfo_full(FILE *f_out, TElliptical *elliptical);
PRECISION VarianceScale_TruncatedPower(int dim, PRECISION lower_bound, PRECISION upper_bound, PRECISION power);
PRECISION VarianceScale_TruncatedGaussian(int dim, PRECISION lower_bound, PRECISION upper_bound, PRECISION s);
+PRECISION VarianceScale_Step(int dim, PRECISION *t, int m);
TElliptical* CreateElliptical_TruncatedGaussian(int dim, TVector mean, TMatrix variance, PRECISION lower_bound, PRECISION upper_bound, PRECISION s);
TElliptical* CreateElliptical_Gaussian(int dim, TVector mean, TMatrix variance);
@@ -93,8 +94,6 @@ TElliptical* CreateElliptical_TruncatedPower(int dim, TVector center, TMatrix sc
PRECISION lower_bound, PRECISION upper_bound, PRECISION power);
TElliptical* CreateElliptical_Step(int dim, TVector center, TMatrix scale, PRECISION *table, int m);
-// Testing routines
-void TestEllipticalVariance(void);
#endif
/********************************************************************************
@@ -128,7 +127,7 @@ This gives that the log of the density of z is given by
-ln(abs(det(S)))-ln(2)-0.5*n*ln(pi)+ln(gamma(n/2))-(n-1)*ln(r)+ln(f(r))
-where r=Inverse(S)*(z-c).
+where r=norm(Inverse(S)*(z-c)).
---------------------------------------------------------------------------------
Note that
@@ -153,7 +152,7 @@ The mapping to the data structure is:
center = c
scale = S
quadratic_form = Inverse(S*S')
- logabsdet = -log(abs(det(S)))
+ logabsdet = -ln(abs(det(S)))
dim = n
---------------------------------------------------------------------------------
--
GitLab