Use of maxima in 'statistical engineering'

FWIW Maxima can do this computation piecewise for 25 vars using the boxcar function and using pw.mac in about < 7 minutes.  pw.mac is a package I wrote and is available on my site. http://mysite.verizon.net/res11w2yb/id2.html

(%i1) (load(pw),showtime:true);
Evaluation took 0.0000 seconds (0.0000 elapsed)
(out1) true
(%i2) convolution(f, g, x) := block([s,p], f:subst(x = p, f), g:subst(x = x - p, g), pwdefint(f * g, p, minf, inf))$
Evaluation took 0.0000 seconds (0.0000 elapsed)
(%i3) define(boxcar(x), 1/(2*a)*(unit_step(x+a)-unit_step(x-a)));
Evaluation took 0.0000 seconds (0.0000 elapsed)
(out3) boxcar(x):=(unit_step(x+a)-unit_step(x-a))/(2*a)
(%i4) dc:convolution(boxcar(x),boxcar(x),x);
Evaluation took 0.0800 seconds (0.0800 elapsed)
(out4) ((x+2*a)*signum(x+2*a)+(x-2*a)*signum(x-2*a)-2*x*signum(x))/(8*a^2)
(%i5) for i : 1 thru 24 do dc:convolution(dc,boxcar(x),x);
Evaluation took 409.7100 seconds (409.7100 elapsed)
(out5) done
(%6) wxdraw2d(
yrange=[0,.4],
color=black,
explicit(dc,x,-5,5)),wxplot_size=[600,400],a=1/2$


Rich

----- Original Message ----- 
From: Vishal Ramnath 
To: maxima at math.utexas.edu 
Sent: Tuesday, August 11, 2009 5:45 AM
Subject: Use of maxima in 'statistical engineering'


Dear all,

Firstly apologies for possibly repeating some of the discussions that have already taken place.

I'm not too sure (I have been on leave on the last week) but I think Robert Dodier mentioned that Maxima has found some application at NIST

http://math.nist.gov/mcsd/

but I can't seem to find the subject thread where he mentioned how Maxima is being used by the mathematical/statistical scientists over there. For those people interested in using Maxima for more statistical work I would recommend having a look at the corresponding British institute

http://www.npl.co.uk/mathematics-scientific-computing/software-support-for-metrology/

which I find to be quite useful and helpful in my work.

Just a brief comment on the "Ideas welcome" subject thread started by Jack Schmidt on the use of Maxima for statistical analysis and a bit on the subject thread on doing convolution integrals.

In my work I have to compute the uncertainties of various mechanical measurements at national laboratory level and I frequently make use of Maxima (symbolic + arbitrary precision) because Excel can oftentimes be a waste of time for more precise/advanced calculations.

In this field of work use is made of some applications of Bayesian statistics and the central limit theorem to approximate the probability distribution of the final underlying measurement (which is defined exactly as a convolution integral via. the Markov formula).

To do the calculations for more advanced and critical measurements one must typically handle 25+ random variables in the convolution integral. I have done a Maxima implementation of the Markov convolution integral that is used in metrology measurement uncertainty work and my experience with the probability density functions that are typically utilized in metrology i.e. Gaussian, rectangular, triangular, U-shaped etc. is that to actually compute the convolution integral one must truncate the pdf's.

Trying to compute the convolution integral where the integration limits run from -\infty to +\infty would frequently result in no solution that Maxima can compute. I don't have access to Mathematica but I am guessing that Mathematica will also have difficulties directly handling the convolution integrals unless the pdf's are truncated (piece-wise continuous is not absolutely essential). 

At this point in time there is also a slight shift amongst some metrology R&D people to use actual pdf's (from experimental as well as simulation data) and not just curve fits of typical pdf's like Gaussian distributions i.e. using "messy" real life pdf's, which means doing a convolution integral to get an answer becomes more difficult. 

To get around this one is then forced to use a Monte Carlo simulation. I have run both programs using open source environments on the same test problem on my office computer using the binary distributions of Maxima (v 5.18.1) for the direct symbolic calculations of the convolution integral and Gnu Octave (v 3.1.50) for a Monte Carlo simulation (using the built in Mersenne twister generator to randomly sample) and my opinion is that at this point in time, at least for metrology work at national lab level, that a MC approach is more feasible and realistic that an attempted direct attack on the convolution integral.

Maxima will work on truncated limits but would for the most part take too long to converge, and will in certain cases fail to find an answer. A Monte Carlo statistical analysis although requiring some 1E5 to 1E6 simulation events will always yield a answer, the only problem is that a lot of computing power is then necessary. 

I have no idea how to use multi-processors like dual core / quad core PC's to get and speed up symbolic computations in Maxima, but the corresponding usage of multi-processors for Monte Carlo work is much more simple since the simulation is embarrassingly parallel. 

Kind regards,




      Vishal Ramnath  
      R&D Metrologist:  Pressure & Vacuum 
     
      Private Bag X34 
      Lynnwood Ridge 
      0040 
      South Africa 
      Tel: +27 12 841 4344 
      Fax: +27 12 841 4458 
     
      e-mail: vramnath at nmisa.org 
      http: www.nmisa.org 



--------------------------------------------------------------------------------


_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/gif
Size: 10664 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20090811/46ebcee2/attachment-0001.gif