Dear Simon,
Thank you for the quick and informative response! I have a few
questions/observations:
1) Thank you for the information regarding the randn() function. I found a
workaround in the original code but now it looks much better and acts as
expected.
2) When I use A = sparse(In - rho*W) I get the following warning from Julia
and a break:
julia> sar_gibbs
LoadError("C:/Users/Don/Desktop/Julia Materials/SAR
Gibbs/sar_gibbs.jl",13,ErrorException("The inverse of a sparse matrix can often
be dense and can cause the computer to run out of memory. If you are sure you
have enough memory, please convert your matrix to a dense matrix."))
Apparently, it does not like it when I invert the A matrix.
3) Can I wrap the code with a "blank" function call like so: function
sar_gibbs()?
4) Are you indicating that I can include the knn_weights_matrix code in the
same file?
Once again, thank you for the response. My goal is to produce a suite of
Bayesian spatial econometric models for public consumption and you help is
greatly appreciated.
Regards,
Don
On Tuesday, July 8, 2014 5:06:14 AM UTC-4, Simon Byrne wrote:
>
> Hi Don,
>
> I think the reason they're not sparse is that in the line
> A = (In - rho*W);
> the matrix In is not sparse: if you replace the line
> In = eye(n);
> by
> In = speye(n);
> the result should then be sparse. However at the moment there doesn't seem
> to be a sparse det method (I just filed issue #7543
> <https://github.com/JuliaLang/julia/issues/7543>): you can get around
> this by using log(det(lufact(A))) instead of log(det(A)).
>
> You might also want to have a quick look at the performance tips
> <http://docs.julialang.org/en/latest/manual/performance-tips/>, in
> particular:
>
> * You should use randn() instead of randn(1): randn() samples a single
> Float64 scalar, randn(1) samples a vector of length 1. On my machine (using
> 0.3 pre-release) your code doesn't work, as this promotes some vectors to
> matrices, which doesn't work with dot. You also then won't need to call
> variables like rho_proposed[1].
> * Try wrapping everything in a function: this makes it easier for the JIT
> compiler to work its magic, as well as making it easier to use profiling
> tools.
> * If you're using dense matrices, you can use logdet(A) instead of
> log(det(A)): this can avoid some underflow and overflow problems
> (unfortunately this doesn't work for sparse lu factorizations yet, see the
> issue mentioned above).
> * Purely from a style point of view, there's no need to keep functions in
> separate files, or end lines with semi-colons.
>
> Simon
>
>
>
> On Tuesday, 8 July 2014 01:48:56 UTC+1, Donald Lacombe wrote:
>>
>> Dear Julia Users,
>>
>> I am currently developing/converting several Bayesian spatial econometric
>> models from MATLAB into Julia. I have successfully coded the spatial
>> autoregressive model (SAR) with diffuse priors in Julia but have a question
>> regarding the use of sparse matrix algorithms. The models use a spatial
>> weight matrix which is usually sparse in nature and this matrix appears in
>> many of the expressions, especially in the random walk Metropolis algorithm
>> used to obtain the marginal distribution for the spatial autocorrelation
>> parameter rho. Here is a code snippet and the complete code is attached:
>>
>> # Draw for rho
>>
>> A = (In - rho*W);
>>
>> denominator = log(det(A))-.5*sigma^-2*dot(A*y-x*beta,A*y-x*beta);
>>
>> accept = 0;
>>
>> rho_proposed = rho + zta*randn(1);
>>
>> while accept == 0
>>
>> if ((rho_proposed[1] > -1) & (rho_proposed[1] < 1));
>>
>> accept = 1;
>>
>> else
>>
>> rho_proposed = rho + zta*randn(1);
>>
>> end
>>
>> end
>>
>> B = (In - rho_proposed[1]*W);
>>
>> numerator = log(det(B))-.5*sigma^-2*dot(B*y-x*beta,B*y-x*beta);
>>
>> u = rand(1);
>>
>> if ((numerator[1] - denominator[1]) > exp(1))
>>
>> pp = 1;
>>
>> else
>>
>> ratio = exp(numerator[1] - denominator[1]);
>>
>> pp = min(1,ratio);
>>
>> end
>>
>> if (u[1] < pp[1])
>>
>> rho = rho_proposed[1];
>>
>> acc = acc + 1;
>>
>> end
>>
>> ar = acc/gibbs;
>>
>>
>> if ar < 0.4
>>
>> zta = zta/1.1;
>>
>> end
>>
>> if ar > 0.6
>>
>> zta = zta*1.1;
>>
>> end
>>
>>
>> A = (In - rho*W);
>>
>>
>> Basically, I'm wondering if there is any way to make the "A" and "B"
>> matrices sparse and possibly make it run faster, especially in the
>> log(det(A)) terms. Currently, 110 draws (with 10 burn) takes approximately 8
>> seconds on my 64 bit, core i7 laptop. The computational speed decreases with
>> the sample size n because the weight matrices are treated as full.
>>
>>
>> Any help would be greatly appreciated and if anyone is interested in running
>> the code, the Distributions and Distance packages must be included and
>> initiated first.
>>
>>
>> Regards,
>>
>> Don
>>
>>
>>
>>
