4
$\begingroup$

I'm doing a bootstrap based on this thread with the following code (file.txt here):

import = Import["file.txt", "List"];
data = Log[10, import];
data = Delete[data, Position[data, RankedMin[data, 1]][[1, 1]]];
data = Delete[data, Position[data, RankedMin[data, 1]][[1, 1]]];
data = Delete[data, Position[data, RankedMax[data, 1]][[1, 1]]];

dist=MixtureDistribution[{p,1-p},{NormalDistribution[μ1, σ1],NormalDistribution[μ2, σ2]}];
param=FindDistributionParameters[data,dist,{{μ1,-0.02},{σ1,0.66},{μ2,1.47},{σ2,0.47},{p,0.21}},ParameterEstimator->{"MaximumLikelihood",Method->"NMaximize"}];

bootstrap=Table[FindDistributionParameters[RandomVariate[dist/.param,Length[data]],dist,{{μ1,-0.02},{σ1,0.66},{μ2,1.47},{σ2,0.47},{p,0.21}},ParameterEstimator->{"MaximumLikelihood",Method->"NMaximize"}],{10}];

param takes 11 sec, I'd expect bootstrap to take approx. 110 sec or something a bit longer (RandomVariate itself takes 0.02 sec, so seems irrelevant for the absolute timing), but actually it takes about 400 sec - looks like $O(n^2)$.

On a 4-core machine, using ParallelTable I can get down to about 100 sec, but that's a side issue - I want to do it 100 or 1000 times. In the latter I expect the code to take $\approx 11\times 1000 sec=3 hrs$, and using ParallelTable go under 1 hr.

In my understanding the above code does the same task (taking 11 sec) multiple times independently of each other, so two tasks should take 22 sec and so on. However, even Table with only two iterations take 36 sec.

Where am I wrong and how to optimize this?

EDIT I checked what happens when I change p to E^w/(1+E^w) (also changing Table to ParallelTable to speed things up):

dist=MixtureDistribution[{E^w/(1+E^w),1/(1+E^w)},{NormalDistribution[μ1, σ1],NormalDistribution[μ2, σ2]}];
param=FindDistributionParameters[data,dist,{{μ1,-0.02},{σ1,0.66},{μ2,1.47},{σ2,0.47},{w,Log[1/3.]}},ParameterEstimator->{"MaximumLikelihood",Method->"NMaximize"}];

bootstrap=ParallelTable[FindDistributionParameters[RandomVariate[dist/.param,Length[data]],dist,{{μ1,-0.02},{σ1,0.66},{μ2,1.47},{σ2,0.47},{w,Log[1/3.]}},ParameterEstimator->{"MaximumLikelihood",Method->"NMaximize"}],{10}];

and it takes 162 sec compared to 115 sec using p as in my initial post.

EEDIT2 After Sjoerd's answer I checked what happens when I set in my initial bootstrap code {p,0.4}: 100 executions take 230 sec, 200 executions - 470 sec. Almost linearly. So it appears, like in the accepted answer, that near the correct values the algorithms needs many iterations to converge and it takes lots of time. Probably setting also other starting values (for $\mu_i$ and $\sigma_i$) would improve the timing even more.

$\endgroup$
  • 1
    $\begingroup$ The variable import on the first line should probably be data, right? $\endgroup$ – Sjoerd C. de Vries Jun 5 '15 at 14:08
  • 1
    $\begingroup$ For a better way to estimate the parameters of a MixtureDistribution have a look at this answer. The current method to determine the weights (p and 1-p) is troublesome and this answer has a better way. Using this simple change makes it already 4 times faster. $\endgroup$ – Sjoerd C. de Vries Jun 5 '15 at 14:15
  • 1
    $\begingroup$ Multiple executions of FindDistributionParameters take about 25% more than based on single execution time, at least in my measurements with four evaluations. My guess is that this is due to garbage collection that takes place on internal data structures that were used by FindDistributionParameters. If I track CPU usage of a single execution I always see the end of the calculation being followed by a large spikesome time after the calculation ends. $\endgroup$ – Sjoerd C. de Vries Jun 5 '15 at 14:36
  • 1
    $\begingroup$ Did you change the parameter to be estimated from p to w? I did, and my times really divide by 4. I don't make these numbers up. And the point is not whether MixtureDistribution knows the parameter range but whether FindDistributionParameters knows this. It has really no way of knowing this; so, this is where things go wrong. $\endgroup$ – Sjoerd C. de Vries Jun 5 '15 at 14:47
  • 1
    $\begingroup$ It appears the execution time is very sensitive to the starting value of w. I used 0.21, you used Log[1/3.]. The latter is 5.5 times slower than the former on my PC. Starting the p version of the estimation on p=0.55 speeds up the single estimation by a factor of 10. $\endgroup$ – Sjoerd C. de Vries Jun 5 '15 at 15:02
4
$\begingroup$

Collecting timings for a range of starting values for the p parameter shows that the timing is critically dependent on this value.

enter image description here

The peak happens to be close to the fitted value of p, 0.215. I assume that the gradient in the neigborhood is so low that the algorithm needs many iterations to converge. The multidimensionality of the situation won't help either.

$\endgroup$
  • $\begingroup$ Interesting observation. It appears that my initial code with {p,0.4} works even faster than using E^w/(1+E^w). Maybe setting other parameters further from the right values will also speed up the computations. $\endgroup$ – corey979 Jun 5 '15 at 16:00
  • 1
    $\begingroup$ Well, that's good to hear. $\endgroup$ – Sjoerd C. de Vries Jun 5 '15 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.