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.

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ The variable import on the first line should probably be data, right? $\endgroup$
    – Sjoerd C. de Vries
    Commented Jun 5, 2015 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
    Commented Jun 5, 2015 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
    Commented Jun 5, 2015 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
    Commented Jun 5, 2015 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
    Commented Jun 5, 2015 at 15:02

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
2
  • $\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
    Commented Jun 5, 2015 at 16:00
  • 1
    $\begingroup$ Well, that's good to hear. $\endgroup$
    – Sjoerd C. de Vries
    Commented Jun 5, 2015 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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