Why yes, that's much better. Very clever way of filling the triangle.
Henry Rich
On 8/22/2012 4:28 PM, Ric Sherlock wrote:
Note that Ewart Shaw has responded to this on the Wiki with a solution
that only generates the required number of random numbers:
http://www.jsoftware.com/jwiki/EwartShaw/RandomSymmetricMatrix
On Wed, Aug 22, 2012 at 5:30 PM, Ric Sherlock <tikk...@gmail.com> wrote:
The other option is not to add them in the first place?
load 'stats/distribs'
(|: + ~: zeroTri ) >: zeroTri rnorm 5 5
0.346799 _1.22161 0.57274 0.556122 _0.329658
_1.22161 0.149955 _1.77435 _1.76668 0.831557
0.57274 _1.77435 0.77674 _0.0690683 _0.967551
0.556122 _1.76668 _0.0690683 0.720588 _0.195658
_0.329658 0.831557 _0.967551 _0.195658 _0.12314
In other words zero the items above the diagonal then transpose and
add to the non-diagonal items.
Where zeroTri is an adverb defined as below:
NB.*zeroTri a Zeros triangular items of matrix determined by verb to left
NB. EG: < zeroTri mat NB. zeros lower-tri items of mat
NB. EG: <: zeroTri mat NB. zeros lower-tri, off-diag items of mat
NB. EG: > zeroTri mat NB. zeros upper tri items of mat
NB. EG: = zeroTri mat NB. zeros all off-diag items of mat
NB. EG: ~: zeroTri mat NB. zeros diag items of mat
zeroTri=: 1 :'([: u/~ i.@#) * ]'
This is somewhat wasteful in that you generate a bunch of random
numbers that you then discard, but if that was an issue you could work
around it.
On Wed, Aug 22, 2012 at 5:07 PM, Owen Marschall
<omarschal...@amherst.edu> wrote:
You read my mind. I was thinking about this same problem tonight while at the
opera, and I couldn't think of any way to only divide the diagonals by sqrt(2)
a second time--without loops, of course. I don't quite yet understand why your
solution works, but I'm sure with enough staring and dictionary help I'll get
it.
Thanks,
Owen
On Aug 21, 2012, at 8:25 PM, Henry Rich wrote:
If you have a matrix a of standard normal deviates, you can make it symmetric
with
asymm =: (+ |:) a
but what is the variance of the items of a?
The variance of values off the principal diagonal will be the sum of the
variance of two independent standard normal deviates. i.e. 2.
To return these values to variance 1 you need to divide by sqrt(2).
But the variance of values ON the principal diagonal will be the sum of two
perfectly correlated random variables, i. e. 4.
So you need to treat the principal diagonal differently. You can reduce its
variance by scaling it differently after the conversion to symmetric, dividing
the diagonal by sqrt(4) and the rest by sqrt(2):
asymmgood =: asymm % %: +: >: e. i. # asymm
(The e. i. bit is a standard idiom for making an identity matrix. Another one
you see around is = i. but I avoid that because I think monad = was wrongly
defined and should be assigned for other purposes)
If I've made a statistical blunder I'm sure someone will tell me.
Henry Rich
On 8/21/2012 3:10 PM, Owen Marschall wrote:
Ah, just what I needed. Thanks!
On Aug 21, 2012, at 1:06 PM, Ric Sherlock wrote:
The primitive ( |: ) is transpose. E.g. :
|: i. 3 4
0 4 8
1 5 9
2 6 10
3 7 11
On Aug 22, 2012 6:55 AM, "Owen Marschall" <omarschal...@amherst.edu> wrote:
Anyone know of an easy way to create a random symmetric matrix (more
specifically, a matrix whose entires are each picked from a standard
Gaussian distribution)? I can start by doing
load 'stats'
R=:normalrand N N
but this is not symmetric, and I don't know of any way to symmetrize it
without thinking in loops, which I'm training myself not to. If I could
somehow take a transpose, that would solve the problem, but I don't know
how to do that either.
Thanks,
Owen
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
