← Back to branch summary

~ubuntu-branches/ubuntu/maverick/openturns/maverick

« back to all changes in this revision

Viewing changes to lib/src/Uncertainty/StatTests/FittingTest.cxx

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2008-11-18 06:32:22 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20081118063222-pa0qncclrerrqkg2
Tags: 0.12.2-1
http://bugs.debian.org/506005
* New upstream release
* Bug fix: "New upstream release available (0.12.2)", thanks to Jerome
  Robert (Closes: #506005).
* Applied patch by J. Robert.
* debian/control: build-depends on libxml2

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/src/Uncertainty/StatTests/FittingTest.cxx
20
20
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
21
21
 *
22
22
 *  @author: $LastChangedBy: dutka $
23
 
 *  @date:   $LastChangedDate: 2008-06-26 13:50:17 +0200 (jeu, 26 jun 2008) $
24
 
 *  Id:      $Id: FittingTest.cxx 862 2008-06-26 11:50:17Z dutka $
 
23
 *  @date:   $LastChangedDate: 2008-10-31 11:52:04 +0100 (ven 31 oct 2008) $
 
24
 *  Id:      $Id: FittingTest.cxx 995 2008-10-31 10:52:04Z dutka $
25
25
 */
26
26
#include <cmath>
27
27
#include <fstream>
31
31
#include "Path.hxx"
32
32
#include "ResourceMap.hxx"
33
33
#include "Log.hxx"
 
34
#include "DistFunc.hxx"
34
35
 
35
36
namespace OpenTURNS
36
37
{
44
45
      typedef Base::Common::Log          Log;
45
46
      typedef Base::Type::NumericalPoint NumericalPoint;
46
47
      typedef Base::Type::Description    Description;
 
48
      typedef Distribution::DistFunc     DistFunc;
47
49
 
48
50
      /* Default constructor, needed by SWIG */
49
51
      FittingTest::FittingTest()
55
57
      FittingTest::Distribution FittingTest::BestModelBIC(const NumericalSample  & sample,
56
58
                                                          const FactoryCollection & factoryCollection)
57
59
      {
58
 
        UnsignedLong size(factoryCollection.getSize());
 
60
        const UnsignedLong size(factoryCollection.getSize());
59
61
        if (size == 0) throw InternalException(HERE) << "Error: no model given";
60
62
        Distribution bestDistribution(factoryCollection[0].buildImplementation(sample));
61
63
        NumericalScalar bestConcordanceMeasure(BIC(sample, bestDistribution, bestDistribution.getParametersNumber()));
62
 
        for (UnsignedLong i = 1; i < size; i ++)
 
64
        for (UnsignedLong i = 1; i < size; ++i)
63
65
          {
64
 
            Distribution distribution(factoryCollection[i].buildImplementation(sample));
65
 
            NumericalScalar concordanceMeasure(BIC(sample, distribution, distribution.getParametersNumber()));
 
66
            const Distribution distribution(factoryCollection[i].buildImplementation(sample));
 
67
            const NumericalScalar concordanceMeasure(BIC(sample, distribution, distribution.getParametersNumber()));
66
68
            if (concordanceMeasure < bestConcordanceMeasure)
67
69
              {
68
70
                bestConcordanceMeasure = concordanceMeasure;
76
78
      FittingTest::Distribution FittingTest::BestModelBIC(const NumericalSample  & sample,
77
79
                                                          const DistributionCollection & distributionCollection)
78
80
      {
79
 
        UnsignedLong size(distributionCollection.getSize());
 
81
        const UnsignedLong size(distributionCollection.getSize());
80
82
        if (size == 0) throw InternalException(HERE) << "Error: no model given";
81
83
        Distribution bestDistribution(distributionCollection[0]);
82
84
        NumericalScalar bestConcordanceMeasure(BIC(sample, bestDistribution));
83
 
        for (UnsignedLong i = 1; i < size; i ++)
 
85
        for (UnsignedLong i = 1; i < size; ++i)
84
86
          {
85
 
            Distribution distribution(distributionCollection[i]);
86
 
            NumericalScalar concordanceMeasure(BIC(sample, distribution));
 
87
            const Distribution distribution(distributionCollection[i]);
 
88
            const NumericalScalar concordanceMeasure(BIC(sample, distribution));
87
89
            if (concordanceMeasure < bestConcordanceMeasure)
88
90
              {
89
91
                bestConcordanceMeasure = concordanceMeasure;
100
102
      FittingTest::Distribution FittingTest::BestModelKolmogorov(const NumericalSample  & sample,
101
103
                                                                 const FactoryCollection & factoryCollection)
102
104
      {
103
 
        UnsignedLong size(factoryCollection.getSize());
 
105
        const UnsignedLong size(factoryCollection.getSize());
104
106
        if (size == 0) throw InternalException(HERE) << "Error: no model given";
105
 
        NumericalScalar fakeLevel(0.5);
 
107
        const NumericalScalar fakeLevel(0.5);
106
108
        Distribution bestDistribution(factoryCollection[0].buildImplementation(sample));
107
109
        NumericalScalar bestPValue(Kolmogorov(sample, bestDistribution, fakeLevel, bestDistribution.getParametersNumber()).getPValue());
108
 
        for (UnsignedLong i = 1; i < size; i ++)
 
110
        for (UnsignedLong i = 1; i < size; ++i)
109
111
          {
110
 
            Distribution distribution(factoryCollection[i].buildImplementation(sample));
111
 
            NumericalScalar pValue(Kolmogorov(sample, distribution, fakeLevel, distribution.getParametersNumber()).getPValue());
 
112
            const Distribution distribution(factoryCollection[i].buildImplementation(sample));
 
113
            const NumericalScalar pValue(Kolmogorov(sample, distribution, fakeLevel, distribution.getParametersNumber()).getPValue());
112
114
            if (pValue > bestPValue)
113
115
              {
114
116
                bestPValue = pValue;
122
124
      FittingTest::Distribution FittingTest::BestModelKolmogorov(const NumericalSample  & sample,
123
125
                                                                 const DistributionCollection & distributionCollection)
124
126
      {
125
 
        UnsignedLong size(distributionCollection.getSize());
 
127
        const UnsignedLong size(distributionCollection.getSize());
126
128
        if (size == 0) throw InternalException(HERE) << "Error: no model given";
127
129
        Distribution bestDistribution(distributionCollection[0]);
128
130
        NumericalScalar bestPValue(Kolmogorov(sample, bestDistribution).getPValue());
129
 
        for (UnsignedLong i = 1; i < size; i ++)
 
131
        for (UnsignedLong i = 1; i < size; ++i)
130
132
          {
131
 
            Distribution distribution(distributionCollection[i]);
132
 
            NumericalScalar pValue(Kolmogorov(sample, distribution).getPValue());
 
133
            const Distribution distribution(distributionCollection[i]);
 
134
            const NumericalScalar pValue(Kolmogorov(sample, distribution).getPValue());
133
135
            if (pValue > bestPValue)
134
136
              {
135
137
                bestPValue = pValue;
144
146
      FittingTest::Distribution FittingTest::BestModelChiSquared(const NumericalSample  & sample,
145
147
                                                                 const FactoryCollection & factoryCollection)
146
148
      {
147
 
        UnsignedLong size(factoryCollection.getSize());
 
149
        const UnsignedLong size(factoryCollection.getSize());
148
150
        if (size == 0) throw InternalException(HERE) << "Error: no model given";
149
 
        NumericalScalar fakeLevel(0.5);
 
151
        const NumericalScalar fakeLevel(0.5);
150
152
        Distribution bestDistribution(factoryCollection[0].buildImplementation(sample));
151
153
        NumericalScalar bestPValue(ChiSquared(sample, bestDistribution, fakeLevel, bestDistribution.getParametersNumber()).getPValue());
152
 
        for (UnsignedLong i = 1; i < size; i ++)
 
154
        for (UnsignedLong i = 1; i < size; ++i)
153
155
          {
154
 
            Distribution distribution(factoryCollection[i].buildImplementation(sample));
155
 
            NumericalScalar pValue(ChiSquared(sample, distribution, fakeLevel, distribution.getParametersNumber()).getPValue());
 
156
            const Distribution distribution(factoryCollection[i].buildImplementation(sample));
 
157
            const NumericalScalar pValue(ChiSquared(sample, distribution, fakeLevel, distribution.getParametersNumber()).getPValue());
156
158
            if (pValue > bestPValue)
157
159
              {
158
160
                bestPValue = pValue;
166
168
      FittingTest::Distribution FittingTest::BestModelChiSquared(const NumericalSample  & sample,
167
169
                                                                 const DistributionCollection & distributionCollection)
168
170
      {
169
 
        UnsignedLong size(distributionCollection.getSize());
 
171
        const UnsignedLong size(distributionCollection.getSize());
170
172
        if (size == 0) throw InternalException(HERE) << "Error: no model given";
171
173
        Distribution bestDistribution(distributionCollection[0]);
172
174
        NumericalScalar bestPValue(ChiSquared(sample, bestDistribution).getPValue());
173
 
        for (UnsignedLong i = 1; i < size; i ++)
 
175
        for (UnsignedLong i = 1; i < size; ++i)
174
176
          {
175
 
            Distribution distribution(distributionCollection[i]);
176
 
            NumericalScalar pValue(ChiSquared(sample, distribution).getPValue());
 
177
            const Distribution distribution(distributionCollection[i]);
 
178
            const NumericalScalar pValue(ChiSquared(sample, distribution).getPValue());
177
179
            if (pValue > bestPValue)
178
180
              {
179
181
                bestPValue = pValue;
190
192
        throw(InvalidArgumentException)
191
193
      {
192
194
        if (sample.getDimension() != distribution.getDimension()) throw InvalidArgumentException(HERE) << "Error: the sample dimension and the distribution dimension must be equal";
193
 
        UnsignedLong size(sample.getSize());
194
 
        UnsignedLong parametersNumber(distribution.getParametersNumber());
 
195
        const UnsignedLong size(sample.getSize());
 
196
        const UnsignedLong parametersNumber(distribution.getParametersNumber());
195
197
        if (parametersNumber < estimatedParameters) throw InvalidArgumentException(HERE) << "Error: the number of estimated parameters cannot exceed the number of parameters of the distribution";
196
198
        NumericalScalar logLikelihood(0.0);
197
 
        for (UnsignedLong i = 0; i < size; i++)
 
199
        for (UnsignedLong i = 0; i < size; ++i)
198
200
          {
199
201
            logLikelihood += log(distribution.computePDF(sample[i]));
200
202
          }
206
208
                                       const Factory & factory)
207
209
        throw(InvalidArgumentException)
208
210
      {
209
 
        Distribution distribution(factory.buildImplementation(sample));
 
211
        const Distribution distribution(factory.buildImplementation(sample));
210
212
        return BIC(sample, distribution, distribution.getParametersNumber());
211
213
      }
212
214
 
219
221
      {
220
222
        if ((level <= 0.0) || (level >= 1.0)) throw InvalidArgumentException(HERE) << "Error: level must be in ]0, 1[, here level=" << level;
221
223
        if (sample.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: Kolmogorov test works only with 1D samples";
222
 
        Distribution distribution(factory.buildImplementation(sample));
 
224
        const Distribution distribution(factory.buildImplementation(sample));
223
225
        if (!distribution.getImplementation()->isContinuous()) throw InvalidArgumentException(HERE) << "Error: Kolmogorov test can be applied only to a continuous distribution";
224
226
        if (distribution.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: Kolmogorov test works only with 1D distribution";
225
227
        return Kolmogorov(sample, distribution, level, distribution.getParametersNumber());
238
240
        if (!distribution.getImplementation()->isContinuous()) throw InvalidArgumentException(HERE) << "Error: Kolmogorov test can be applied only to a continuous distribution";
239
241
        if (distribution.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: Kolmogorov test works only with 1D distribution";
240
242
        if (estimatedParameters > 0) Log::Info("Warning: using Kolmogorov test for a distribution with estimated parameters will result in an overestimated pValue");
241
 
        return RunRTest(sample, distribution, level, estimatedParameters, "Kolmogorov");
 
243
        const NumericalSample sortedSample(sample.sort(0));
 
244
        const UnsignedLong size(sample.getSize());
 
245
        NumericalScalar value(0.0);
 
246
        for (UnsignedLong i = 0; i < size; ++i)
 
247
          {
 
248
            const NumericalScalar cdfValue(distribution.computeCDF(sortedSample[i]));
 
249
            value = std::max(value, std::max(fabs(NumericalScalar(i) / size - cdfValue), fabs(cdfValue - NumericalScalar(i + 1) / size)));
 
250
          }
 
251
        const NumericalScalar pValue(DistFunc::pKolmogorov(size, value, true));
 
252
        TestResult result(OSS() << "Kolmogorov" << distribution.getClassName(), (pValue > 1.0 - level), pValue, 1.0 - level);
 
253
        return result;
242
254
      }
243
255
 
244
256
      /* Chi-squared test */
249
261
      {
250
262
        if ((level <= 0.0) || (level >= 1.0)) throw InvalidArgumentException(HERE) << "Error: level must be in ]0, 1[, here level=" << level;
251
263
        if (sample.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: ChiSquared test works only with 1D samples";
252
 
        Distribution distribution(factory.buildImplementation(sample));
 
264
        const Distribution distribution(factory.buildImplementation(sample));
253
265
        if (distribution.getImplementation()->isContinuous()) throw InvalidArgumentException(HERE) << "Error: Chi-squared test cannot be applied to a continuous distribution";
254
266
        if (distribution.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: ChiSquared test works only with 1D distribution";
255
267
        return ChiSquared(sample, distribution, level, distribution.getParametersNumber());
278
290
                                                    const String & testName)
279
291
        throw(InternalException)
280
292
      {
281
 
        String dataFileName(sample.storeToTemporaryFile());
282
 
        String resultFileName(Path::BuildTemporaryFileName("RResult.txt.XXXXXX"));
283
 
        String commandFileName(Path::BuildTemporaryFileName("RCmd.R.XXXXXX"));
 
293
        const String dataFileName(sample.storeToTemporaryFile());
 
294
        const String resultFileName(Path::BuildTemporaryFileName("RResult.txt.XXXXXX"));
 
295
        const String commandFileName(Path::BuildTemporaryFileName("RCmd.R.XXXXXX"));
284
296
        std::ofstream cmdFile(commandFileName.c_str(), std::ios::out);
285
297
        // Fill-in the command file
286
298
        cmdFile << "library(rotRPackage)" << std::endl;
289
301
        cmdFile << "sample <- data.matrix(read.table(\"" << dataFileName << "\"))" << std::endl;
290
302
        cmdFile << "res <- computeTest" << testName << distribution.getImplementation()->getClassName();
291
303
        cmdFile << "(sample, ";
292
 
        NumericalPoint parameters(distribution.getParametersCollection()[0]);
293
 
        UnsignedLong parametersNumber(parameters.getDimension());
294
 
        for (UnsignedLong i = 0; i < parametersNumber; i++)
 
304
        const NumericalPoint parameters(distribution.getParametersCollection()[0]);
 
305
        const UnsignedLong parametersNumber(parameters.getDimension());
 
306
        for (UnsignedLong i = 0; i < parametersNumber; ++i)
295
307
          {
296
308
            cmdFile << parameters[i] << ", ";
297
309
          }
303
315
 
304
316
        OSS systemCommand;
305
317
        systemCommand << ResourceMap::GetInstance().get("R-executable-command") << " --no-save --silent < " << commandFileName << " 2>&1 > /dev/null";
306
 
        int returnCode(system(String(systemCommand).c_str()));
 
318
        const int returnCode(system(String(systemCommand).c_str()));
307
319
        if (returnCode != 0) throw InternalException(HERE) << "Error: unable to execute the system command " << String(systemCommand) << " returned code is " << returnCode;
308
320
        // Parse result file
309
321
        std::ifstream resultFile(resultFileName.c_str(), std::ios::in);