Community Tip - If community subscription notifications are filling up your inbox you can set up a daily digest and get all your notifications in a single email. X
I am still using MC 11.
I was able to make a random square matrix called A of n dimension, with 1<n<1000 (infinite would have been better but my pc horse power is limited). A needed to be positive, with Max Eigenvalue lower than 1.
I used the following combination of formulas, which probably is not the optimum, but it works.
1): Z(z,z)=:rnd(1/n) # so to have numbers that fulfill the condition MaxEval<1
2) A=: matrix(n,n,Z)
I would like to refine it, and I ask kindly advice.
The constrains that I would like to add in generating A are the following.
a) Still a nxn square Matrix with Max Eigenvalue lower than 1.
b) But A should be semi-positive, with random coefficients, including zeros and not just positive as above.
c) However A should be indecomposable (called also irreducible), i.e. A^x (with x a positive number) should turn a positive matrix.
Thanks for any suggestion.
Solved! Go to Solution.
With the rnd() random number generator you get random numbers from a uniform sitribution. The chance of hitting a specific number (such as 0) is virtually zero; exactly it is 1 over the number of possible random values, since Mathcad's range of real number is pretty large, you chance of hitting 0 therefor is pretty low: approximatly 0. If you want a specific probability of occurrence of any number, you have to restrict the number of possible random values. I admit that I restricted it fiercely. You may want to use:
Z(z,z)=: floor(rnd(n^2))/n^3
or any higher power
I considered another option to include 0, which is
Z(z,z)=: ((rnd(2*n)-n)^2)/n^2
This brings 0 to the middle of the range, but does not increase the probability of hitting it exactly. On first run with n=4, I got 3.33e-6 as smallest value.
Yet another possibility is to additionally declare all values within a certain distance from 0 as 0 itself:
Z(z,z)=:
|r<-rnd(n)/n^2
|return if(r<0.0001,0,r)
But this ruins your uniform distribution.
Another construct could be:
Z(z,z)=: if(rnd(n)<0.0001,0,rnd(r))/n^2
Which could result in non-zero numbers below 0.0001/n^2, because if on first call to rnd() the result is >= 0.0001, a new random number is generated by the second call to rnd(), but I don't exactly know what it does to your (uniform) distribution.
Success!
Luc
Many Thanks Luc!
Seems simple and working. Yet I noticed that it generates only n numbers, zeros included. So for a matrix 3x3 there are only 3 numbers generated, which are put in the 9 cells. And this may be a problem, since each cell does not have a random semipositive number (I should probably mentioned that in my constraints). I fear that matrix A could become singular, which is not what I want.
With the rnd() random number generator you get random numbers from a uniform sitribution. The chance of hitting a specific number (such as 0) is virtually zero; exactly it is 1 over the number of possible random values, since Mathcad's range of real number is pretty large, you chance of hitting 0 therefor is pretty low: approximatly 0. If you want a specific probability of occurrence of any number, you have to restrict the number of possible random values. I admit that I restricted it fiercely. You may want to use:
Z(z,z)=: floor(rnd(n^2))/n^3
or any higher power
I considered another option to include 0, which is
Z(z,z)=: ((rnd(2*n)-n)^2)/n^2
This brings 0 to the middle of the range, but does not increase the probability of hitting it exactly. On first run with n=4, I got 3.33e-6 as smallest value.
Yet another possibility is to additionally declare all values within a certain distance from 0 as 0 itself:
Z(z,z)=:
|r<-rnd(n)/n^2
|return if(r<0.0001,0,r)
But this ruins your uniform distribution.
Another construct could be:
Z(z,z)=: if(rnd(n)<0.0001,0,rnd(r))/n^2
Which could result in non-zero numbers below 0.0001/n^2, because if on first call to rnd() the result is >= 0.0001, a new random number is generated by the second call to rnd(), but I don't exactly know what it does to your (uniform) distribution.
Success!
Luc