Skip to content

Commit 9411ee2

Browse files
jennyfarmerjennyfarmer
authored
Add files via upload
1 parent 38593c8 commit 9411ee2

30 files changed

+3957
-0
lines changed

‎ChebyShev.cpp‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
/*
2+
* PDF Estimator: A non-parametric probability density estimation tool based on maximum entropy
3+
* File: ChebyShev.cpp
4+
* Copyright (C) 2018
5+
* Jenny Farmer [email protected]
6+
* Donald Jacobs [email protected]
7+
*
8+
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published
9+
* by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in
10+
* the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11+
* PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with
12+
* this program. If not, see <http://www.gnu.org/licenses/>.
13+
*/
14+
15+
#include "ChebyShev.h"
16+
17+
ChebyShev::ChebyShev() {
18+
19+
}
20+
21+
ChebyShev::ChebyShev(const ChebyShev& orig) {
22+
}
23+
24+
ChebyShev::~ChebyShev() {
25+
}
26+
27+
void ChebyShev::initialize(double dzLocal[], int sizeLocal) {
28+
this->size = sizeLocal;
29+
this->dz = dzLocal;
30+
vector <double> zeroT;
31+
vector <double> oneT;
32+
vector <double> twoT;
33+
34+
for (int z = 0; z < size; z++) {
35+
double x = -1 + dz[z]*2;
36+
zeroT.push_back(1);
37+
oneT.push_back(x);
38+
twoT.push_back(2*x*oneT[z] - 1);
39+
}
40+
termsT.push_back(zeroT);
41+
termsT.push_back(oneT);
42+
termsT.push_back(twoT);
43+
}
44+
45+
double* ChebyShev::getTerms(unsigned mode) {
46+
if (termsT.size() <= mode) {
47+
vector <double> test = addMode(mode);
48+
return &test[0];
49+
} else {
50+
return &termsT[mode][0];
51+
}
52+
}
53+
54+
vector < vector < double > > ChebyShev::getAllTerms(unsigned mode) {
55+
for (unsigned i = 0; i < mode; i++) {
56+
if (termsT.size() <= i) {
57+
addMode(i);
58+
}
59+
}
60+
return termsT;
61+
}
62+
63+
64+
vector <double> ChebyShev::addMode(int mode) {
65+
66+
vector <double> T = termsT.at(mode-1);
67+
vector <double> Tprev = termsT.at(mode-2);
68+
vector <double> Tnext;
69+
double x = 0;
70+
71+
for (int z = 0; z < size; z++) {
72+
x = -1 + 2*dz[z];
73+
Tnext.push_back(2*x*T[z] - Tprev[z]);
74+
}
75+
termsT.push_back(Tnext);
76+
return Tnext;
77+
}
78+

‎ChebyShev.h‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
/*
2+
* PDF Estimator: A non-parametric probability density estimation tool based on maximum entropy
3+
* File: ChebyShev.h
4+
* Copyright (C) 2018
5+
* Jenny Farmer [email protected]
6+
* Donald Jacobs [email protected]
7+
*
8+
* This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published
9+
* by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in
10+
* the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11+
* PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with
12+
* this program. If not, see <http://www.gnu.org/licenses/>.
13+
*/
14+
15+
#ifndef CHEBYSHEV_HPP
16+
#define CHEBYSHEV_HPP
17+
18+
#include <vector>
19+
20+
using namespace std;
21+
class ChebyShev {
22+
public:
23+
ChebyShev();
24+
ChebyShev(const ChebyShev& orig);
25+
virtual ~ChebyShev();
26+
void initialize(double dzLocal[], int sizeLocal);
27+
double * getTerms(unsigned mode);
28+
vector < vector < double > > getAllTerms(unsigned mode);
29+
30+
private:
31+
int size;
32+
double * dz;
33+
vector < vector<double> > termsT;
34+
35+
vector <double> addMode(int mode);
36+
37+
};
38+
39+
#endif /* CHEBYSHEV_HPP */
40+

‎CompilePDF.m‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
mex -O 'EstimatePDF.cpp' 'OutputControl.cpp' 'callPDF.cpp' 'WriteResults.cpp' 'Score.cpp' 'ScoreLL.cpp' 'ScoreQZ.cpp' 'MinimizeScore.cpp' 'InputParameters.cpp' 'InputData.cpp' 'ChebyShev.cpp'

‎EstimatePDF.cpp‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
/*
2+
3+
* File: PDFMainMatlab.cpp
4+
* Author: jenny
5+
*
6+
* Created on December 3, 2018, 8:32 AM
7+
*/
8+
9+
10+
#include "EstimatePDF.h"
11+
#include "callPDF.h"
12+
13+
14+
15+
void MexFunction::operator()(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs) {
16+
17+
int nParameters = inputs.size();
18+
matlabPtr = getEngine();
19+
out.debug = true;
20+
21+
matlab::data::ArrayFactory factory;
22+
if (nParameters < 1) {
23+
out.displayError(matlabPtr, "Data sample required");
24+
}
25+
if (nParameters > 2) {
26+
out.displayError(matlabPtr, "Only two input arguements allowed: data sample and parameter structure");
27+
}
28+
29+
TypedArray<double> doubleArray = std::move(inputs[0]);
30+
int sampleSize = doubleArray.getNumberOfElements();
31+
if (sampleSize < 3) {
32+
out.displayError(matlabPtr, "Must have at least three data points in sample");
33+
}
34+
double sampleData[sampleSize];
35+
int i = 0;
36+
for (auto& elem : doubleArray) {
37+
sampleData[i++] = elem;
38+
}
39+
40+
callPDF *pd = new callPDF(matlabPtr);
41+
42+
if (nParameters > 1) {
43+
StructArray const matlabStructArray = inputs[1];
44+
auto fields = matlabStructArray.getFieldNames();
45+
std::vector<std::string> fieldNames(fields.begin(), fields.end());
46+
int count = 0;
47+
Array structField;
48+
for (std::vector<string>::iterator iter = fieldNames.begin(); iter != fieldNames.end(); ++iter) {
49+
string field = *iter;
50+
if (strcmp(field.c_str(), "SURDtarget") == 0) {
51+
structField = matlabStructArray[0][fieldNames[count++]];
52+
pd->setSURDtarget(structField[0]);
53+
} else if (strcmp(field.c_str(), "SURDmin") == 0) {
54+
structField = matlabStructArray[0][fieldNames[count++]];
55+
pd->setSURDmin(structField[0]);
56+
} else if (strcmp(field.c_str(), "SURDmax") == 0) {
57+
structField = matlabStructArray[0][fieldNames[count++]];
58+
pd->setSURDmax(structField[0]);
59+
} else if (strcmp(field.c_str(), "scoreType") == 0) {
60+
structField = matlabStructArray[0][fieldNames[count++]];
61+
pd->setScoreType((CharArray(structField)).toAscii());
62+
} else if (strcmp(field.c_str(), "LagrangeMin") == 0) {
63+
structField = matlabStructArray[0][fieldNames[count++]];
64+
pd->setLagrangeMin(structField[0]);
65+
} else if (strcmp(field.c_str(), "LagrangeMax") == 0) {
66+
structField = matlabStructArray[0][fieldNames[count++]];
67+
pd->setLagrangeMax(structField[0]);
68+
} else if (strcmp(field.c_str(), "integrationPoints") == 0) {
69+
structField = matlabStructArray[0][fieldNames[count++]];
70+
pd->setPoints(structField[0]);
71+
} else if (strcmp(field.c_str(), "lowBound") == 0) {
72+
structField = matlabStructArray[0][fieldNames[count++]];
73+
pd->setLow(structField[0]);
74+
} else if (strcmp(field.c_str(), "highBound") == 0) {
75+
structField = matlabStructArray[0][fieldNames[count++]];
76+
pd->setHigh(structField[0]);
77+
} else if (strcmp(field.c_str(), "outlierCutoff") == 0) {
78+
structField = matlabStructArray[0][fieldNames[count++]];
79+
pd->setOutlierCutoff(structField[0]);
80+
} else if (strcmp(field.c_str(), "partition") == 0) {
81+
structField = matlabStructArray[0][fieldNames[count++]];
82+
pd->setPartitionSize(structField[0]);
83+
} else if (strcmp(field.c_str(), "debug") == 0) {
84+
structField = matlabStructArray[0][fieldNames[count++]];
85+
pd->setDebug(structField[0]);
86+
} else if (strcmp(field.c_str(), "minVariance") == 0) {
87+
structField = matlabStructArray[0][fieldNames[count++]];
88+
pd->setVariance(structField[0]);
89+
} else if (strcmp(field.c_str(), "adaptiveDx") == 0) {
90+
structField = matlabStructArray[0][fieldNames[count++]];
91+
pd->setAdaptiveDx(structField[0]);
92+
} else {
93+
string unknown = "Unknown parameter: " + field;
94+
out.displayError(matlabPtr, unknown);
95+
}
96+
}
97+
}
98+
pd->makeCall(sampleSize, sampleData);
99+
100+
vector <int> failed;
101+
failed.push_back(pd->solutionFailed);
102+
ArrayDimensions sz = {failed.size(), 1, 1};
103+
TypedArray<int> MsolutionFailed = factory.createArray(sz, failed.begin(), failed.end());
104+
outputs[0] = MsolutionFailed;
105+
106+
if (!pd->solutionFailed) {
107+
sz = {pd->Vx.size(), 1, 1};
108+
TypedArray<double> Mx = factory.createArray(sz, pd->Vx.begin(), pd->Vx.end());
109+
outputs[1] = Mx;
110+
sz = {pd->Vpdf.size(), 1, 1};
111+
TypedArray<double> Mpdf = factory.createArray(sz, pd->Vpdf.begin(), pd->Vpdf.end());
112+
outputs[2] = Mpdf;
113+
sz = {pd->Vcdf.size(), 1, 1};
114+
TypedArray<double> Mcdf = factory.createArray(sz, pd->Vcdf.begin(), pd->Vcdf.end());
115+
outputs[3] = Mcdf;
116+
sz = {pd->Vsqr.size(), 1, 1};
117+
TypedArray<double> Msqr = factory.createArray(sz, pd->Vsqr.begin(), pd->Vsqr.end());
118+
outputs[4] = Msqr;
119+
sz = {pd->Vlagrange.size(), 1, 1};
120+
TypedArray<double> Mlagrange = factory.createArray(sz, pd->Vlagrange.begin(), pd->Vlagrange.end());
121+
outputs[5] = Mlagrange;
122+
123+
vector <double> threshold;
124+
threshold.push_back(pd->thresholdScore);
125+
sz = {threshold.size(), 1, 1};
126+
TypedArray<double> Mthreshold = factory.createArray(sz, threshold.begin(), threshold.end());
127+
outputs[6] = Mthreshold;
128+
129+
vector <double> confidence;
130+
confidence.push_back(pd->confidenceScore);
131+
sz = {confidence.size(), 1, 1};
132+
TypedArray<double> Mconfidence = factory.createArray(sz, confidence.begin(), confidence.end());
133+
outputs[7] = Mconfidence;
134+
135+
sz = {pd->Vr.size(), 1, 1};
136+
TypedArray<double> Mr = factory.createArray(sz, pd->Vr.begin(), pd->Vr.end());
137+
outputs[8] = Mr;
138+
}
139+
delete pd;
140+
141+
}

‎EstimatePDF.h‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
/*
2+
* File: EstimatePDF.h
3+
* Author: jenny
4+
*
5+
* Created on February 2, 2019, 2:00 PM
6+
*/
7+
8+
#ifndef ESTIMATEPDF_H
9+
#define ESTIMATEPDF_H
10+
11+
#include "mexAdapter.hpp"
12+
#include "OutputControl.h"
13+
14+
using namespace matlab::data;
15+
16+
class MexFunction : public matlab::mex::Function {
17+
18+
private:
19+
std::shared_ptr<matlab::engine::MATLABEngine> matlabPtr;
20+
OutputControl out;
21+
22+
public:
23+
void operator()(matlab::mex::ArgumentList outputs, matlab::mex::ArgumentList inputs);
24+
};
25+
26+
#endif

‎EstimatePDF.m‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
% EstimatePDF computes a probability density estimate for a one-dimensional
2+
% data sample based on a nonparametric maximum entropy method. For additional
3+
% usage information see readme.txt file with this installation.
4+
% For conceptual details of the method and algorithm, see:
5+
%
6+
% Farmer, Jenny and Donald Jacobs (2018).
7+
% "High throughput nonparametric probability density estimation." PLoS One 13(5): e0196937.
8+
%
9+
%
10+
% [FAILED, XI, F, CDF, SQR, LAGRANGE, SCORE, CONFIDENCE, SURD] = EstimatePDF(X)
11+
% Computes the density estimate of data in sample X with default settings.
12+
% Outputs returned are as follows:
13+
%
14+
% FAILED non-zero if solution not found
15+
% XI relative spacing along x-axis for density estimate
16+
% F probability density function for estimate (PDF)
17+
% CDF cummulative density function for estimate
18+
% SQR scaled quantile residual for sample data
19+
% LAGRANGE Lagrange multipliers estimated to construct PDF
20+
% SCORE Value of final score for estimate
21+
% CONFIDENCE Confidence threshold for estimate
22+
% SURD sampled uniform random data
23+
%
24+
% [...] = EstimatePDF(X, PARM) Computes the density estimate of data in
25+
% sample X with parameters in a structure defined in PARM. Parameter options
26+
% are listed below.
27+
%
28+
% Parameter name Default value Description
29+
%
30+
% lowBound calculated set fixed lower bound
31+
% highBound calculated set fixed upper bound
32+
% integrationPoints max(1500, 200/n + 200) data resolution; n=sample size
33+
% LagrangeMin 1 minimum number of weighted expansions
34+
% LagrangeMax 200 maximum number of weighted expansions
35+
% SURDtarget 40 target confidence threshold
36+
% SURDmin 5 minimum conficence threshold accepted
37+
% SURDmax 100 maximum conficence threshold accepted
38+
% scoreType ‘QZ' Change to ‘LL’ for log likelihood scoring
39+
% debug false detailed output to console
40+
% partition 1025 Initial data partition for scoring, set to zero for no partitioning
41+
% outlierCutoff 7 Set to zero to disable outlier detection
42+
% adaptiveDx true Set to false for uniformly spaced numerical solution
43+
%
44+
%
45+
% Example: Plot estimated density for a Normal distribution using default
46+
% settings:
47+
% [failed, y, pdf] = EstimatePDF(randn(1000, 1));
48+
% if ~failed
49+
% plot(y, pdf);
50+
% end
51+
%
52+
% Example: Plot estimated cummulative density for a uniform distribution defined on
53+
% the interval (0 1):
54+
% parms.lowBound = 0;
55+
% parms.highBound = 1;
56+
% [failed, y, pdf, cdf] = EstimatePDF(rand(1000, 1), parms);
57+
% if ~failed
58+
% plot(y, cdf);
59+
% end
60+
%
61+
% Example: Plot estimated density for an exponential distribution defined on
62+
% the interval (0 inf), requiring at least 2 Langrange multipliers:
63+
% parms.lowBound = 0;
64+
% parms.LagrangeMin = 2;
65+
% [failed, y, pdf] = EstimatePDF(exprnd(1, 1000, 1), parms);
66+
% if ~failed
67+
% plot(y, pdf);
68+
% end
69+
70+
71+
72+
73+

‎FigureSettings.m‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
function FigureSettings()
2+
%FigureSettings creates publication-quality default settings for all
3+
%subsequent plots in current MATLAB session.
4+
set(0,'DefaultFigureColor','white')
5+
fig.InvertHardcopy = 'off';
6+
width = 6; % Width in inches
7+
height = 4; % Height in inches
8+
alw = 1.5; % AxesLineWidth
9+
fsz = 14; % Fontsize
10+
lw = 1.5; % LineWidth
11+
msz = 8; % MarkerSize
12+
set(0,'defaultAxesFontSize',fsz);
13+
set(0,'defaultLineLineWidth',lw);
14+
set(0,'defaultLineMarkerSize',msz);
15+
set(0,'defaultAxesLineWidth',alw);
16+
defpos = get(0,'defaultFigurePosition');
17+
set(0,'defaultFigurePosition', [defpos(1) defpos(2) width*100, height*100]);
18+
set(0,'defaultFigurePosition', [400, 50, width*100, height*110]);
19+
end

0 commit comments

Comments
 (0)