On Tue, 19 Jan 2016, Riccardo (Jack) Lucchetti wrote:
On Sat, 16 Jan 2016, Allin Cottrell wrote:
> My test was admittedly kinda silly, in that there's not really any reason
> to delegate to a "foreign" program stuff that gretl handles well natively.
> One would be more likely to get Julia to do MCMC or the like.
>
> That said, among the various "foreign" programs on which I've tried the
> notorious Longley exercise, only numpy comes close to gretl for numerical
> precision. However, Anders makes a fair point in saying that the
> statistical error (and I would add, data error) swamps the numerical error
> for this sort of linear problem.
This is problaby more akin to a real-life case:
<hansl>
matrix X = mnormal(6,6)
mwrite(X, "X.mat")
foreign language=julia
function ReadGretlMat(fname)
moo = readdlm(fname);
rows = moo[1,1];
cols = moo[1,2];
v = Array(Float64,rows,cols);
[v[i] = convert(Float64,moo'[i+cols]) for i=1:cols*rows]
return(v');
end
function WriteGretlMat(A, fname)
r,c = size(A)
v = Array(Any,(r+1)*c);
v[1] = r;
v[2] = c;
[v[i] = "" for i=3:c];
[v[i+c] = A'[i] for i=1:r*c];
writedlm(fname,reshape(v,c,r+1)');
end
X = ReadGretlMat("X.mat");
WriteGretlMat(inv(X), "invX.mat");
end foreign
iX = mread("invX.mat")
check = X * iX
print check
</hansl>
(with apologies for my poor and non-idiomatic use of Julia --- but I'll
improve, promise).
Semicolons to end statements are not idiomatic in Julia (unless you
want more than one statement per line). Other than that, looks fine.
For comparison, here's the "gretl_io.jl" that gretl writes for
sourcing by julia under "foreign":
<julia>
gretl_dotdir = <whatever>
function gretl_export(M, fname, autodot=1)
r,c = size(M)
if autodot != 0
fname = gretl_dotdir * fname
end
f = open(fname, "w")
@printf(f, "%d\t%d\n", r, c)
for i = 1:r
for j = 1:c
@printf(f, "%.18e ", M[i,j])
end
@printf(f, "\n")
end
close(f)
end
function gretl_loadmat(fname, autodot=1)
if autodot != 0
fname = gretl_dotdir * fname
end
M = readdlm(fname, skipstart=1)
end
</julia>
My gretl_export() is a bit verbose, but gretl_loadmat() works fine
just by discarding the first row of the gretl .mat file.
Allin