I am trying to implement a simple statistical model for the first time in Mathematica.

My model is made my one random variable which is the sum of background and signal. Background is an exponential random variable, while the signal is a gaussian peak. The only free parameters in the model are the number of signal s and background b events and the position of the gaussian mH.

signalModel = TruncatedDistribution[{100., 150.}, NormalDistribution[mH, 1]];
bkgModel = TruncatedDistribution[{100., 150.}, ExponentialDistribution[1/100.]];
model = MixtureDistribution[{s,b}, {signalModel, bkgModel}];

signalExample = 50.; backgroundExample = 5000.;
data = RandomVariate[Evaluate[model /. {s->signalExample, b->backgroundExample, mH->125.}], 5050];

EstimatedDistribution[data, model, {{s, 0}, {b, 4000}, {mH, 125}}] // Timing

This output:


This takes 10 s and the solution is not so good (26 instead of 50, 3760 instead of 5000)

As a comparison I have implemented the same model with RooFit. The only difference is that RooFit creates an extended likelihood. It takes 0.2 s. RooFit is quite smart since it caches partial results when computing the likelihood.

Is it possibile to speedup the Mathematica version? Am I doing something wrong?

In addition I would like to find the solution only in a certain region (rectangular cuts), but I haven't found a way to implement constraints in EstimatedDistribution.


Your Answer

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

Browse other questions tagged or ask your own question.