Question: Program using Statistics, Sample


Please help my with my program. First, I'll post the relevant code, then the questions:


with(Statistics); with(LinearAlgebra);

#k is a Vector everywhere in this program
ModeFreq := (k, m)->Norm(<Norm(k, 2), m>, 2);
OscillatorAmplitude :=(n, k, m, x) -> (ModeFreq(k, m)/Pi)^(1/4)*exp(-1/2*x^2*ModeFreq(k, m))*HermiteH(n, x*sqrt(ModeFreq(k, m)))/sqrt(2^n*factorial(n));
DeltaApproxAmplitude :=(x0, N, k, m, x) -> evalf(expand((sum(OscillatorAmplitude(n, k, m, x)*OscillatorAmplitude(n, k, m, x0), n = 1 .. N))/sqrt(sum(OscillatorAmplitude(n, k, m, x0)^2, n = 1 .. N))));
DeltaApproxPDF := (x0, N, k, m, x) ->DeltaApproxAmplitude(x0, N, k, m, x)^2;
DeltaApproxDist := (x0, N, k, m) -> Distribution(PDF = (x-> DeltaApproxPDF(x0, N, k, m, x)));

#Approximation of a delta distribution peaked at x0 using N terms.
DeltaApproxRV := (x0, N, k, m) -> RandomVariable(DeltaApproxDist(x0, N, k, m));

fudge:= proc(coefficientmatrix, N) local i, j, newmatrix;
newmatrix := Matrix(Dimension(coefficientmatrix));
for i from 1 to RowDimension(coefficientmatrix) do
for j from 1 to ColumnDimension(coefficientmatrix) do
newmatrix[i, j] := Sample(DeltaApproxRV(coefficientmatrix[i, j], N, <i,j>, 0), 1)[1];
end do end do; return newmatrix;
end proc;


Here are my problems:

1) Sampling from DeltaApproxRV is VERY slow. I need to speed it up as much as possible. For example,
time(Sample(DeltaApproxRV(1, 10, <1>, 0), 1)) gives 51.634, and increasing N makes it worse. The fudge procedure on a 15x15 matrix, with N=10 can take 4 hours or more.

2) fudge just doesn't work at all with high enough N. For example, fudge(Matrix(15),20) just gives the following after a few hours:

Error, (in Statistics:-Sample) FAIL .. 2.029919585 is an invalid range

I've also received, on another occasion:

Error, (in Statistics:-Sample) FAIL .. FAIL is an invalid range

I can't seem to replicate this error with individual Sample commands.

Any help will be very appreciated.

Please Wait...