1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
/*=========================================================================
Program: Visualization Toolkit
Module: vtkPCAStatisticsGnuR.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkPCAStatisticsGnuR.h"
#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkRInterface.h"
#include "vtkObjectFactory.h"
#include <sstream>
vtkStandardNewMacro(vtkPCAStatisticsGnuR);
// ----------------------------------------------------------------------
vtkPCAStatisticsGnuR::vtkPCAStatisticsGnuR()
{
}
// ----------------------------------------------------------------------
vtkPCAStatisticsGnuR::~vtkPCAStatisticsGnuR()
{
}
// ----------------------------------------------------------------------
void vtkPCAStatisticsGnuR::PrintSelf( ostream& os, vtkIndent indent )
{
this->Superclass::PrintSelf( os, indent );
}
// Use R to obtain the p-values for the Chi square distribution with 2 DOFs
vtkDoubleArray* vtkPCAStatisticsGnuR::CalculatePValues(vtkIdTypeArray* dimCol,
vtkDoubleArray* statCol)
{
// Prepare VTK - R interface
vtkRInterface* ri = vtkRInterface::New();
// Use the calculated Jarque-Bera-Srivastava statistics as input to the Chi square function
ri->AssignVTKDataArrayToRVariable( statCol, "jbs" );
ri->AssignVTKDataArrayToRVariable( dimCol, "d" );
// Calculate the p-values (p+1 degrees of freedom)
// Now prepare R script and calculate the p-values (in a single R script evaluation for efficiency)
std::ostringstream rs;
rs << "p<-c();"
<< "for(i in 1:"
<< dimCol->GetNumberOfTuples()
<< "){"
<< "p<-c(p,1-pchisq(jbs[i],d[i]));"
<< "}";
ri->EvalRscript( rs.str().c_str() );
// Retrieve the p-values
vtkDoubleArray* testCol = vtkArrayDownCast<vtkDoubleArray>( ri->AssignRVariableToVTKDataArray( "p" ) );
if ( ! testCol || testCol->GetNumberOfTuples() != statCol->GetNumberOfTuples() )
{
vtkWarningMacro( "Something went wrong with the R calculations. Reported p-values will be invalid." );
testCol = this->Superclass::CalculatePValues( dimCol, statCol );
}
else
{
// increment ref count on testCol so its not cleaned up when the R interface goes away
testCol->Register(NULL);
}
// Clean up
ri->Delete();
return testCol;
}
|
