Thanks for your replies. I agree that my code right now is very clunky (and far from efficient) I just wanted to try it out because the R implementation was taking some time.
The basic idea of my code is to implement a linear "unmixing" algorithm to find mixing coefficients (alpha) and their variance sigma2r from multivariate responses with the constraints that the coefficients be positive and sum up to 1. In a 2 component case the coefficients are not independent (i.e. if you know one of them, say alpha1 the other is by default 1 - alpha1. The idea behind the first chunk of code was to choose at random the first coefficient to be estimated each iteration, and at the end of the loop the other is obtained as the difference from 1. I believe I can think of a way to make it more efficient, but I was having trouble with my current implementation. The algorithm is tuned up in a way that the acceptance of each mixing coefficient is around 30-40% for each of them (which is what I've been getting in R on a 2 component simulated data). However, currently (in my Julia implementation) I'm getting high acceptances for either one of them and very low for the other one (85%, 5%), and I wasn't sure what part of my code was not right (besides not being optimal). Cheers, Napo On Sunday, November 30, 2014 7:55:43 AM UTC-6, Johan Sigfrids wrote: > > I think you need to worry about correctness before you worry about speed, > because as it is this code makes no sense to me. For example, at the top of > the code you define: > > R = size(X, 2) > > This means R is an integer. Then, inside your loop, you do: > > tmp = randperm(R) > > Doing randperm() on an integer is the same as doing randperm() on a > vector from 1 to the integer, so now tmp is a vector of the integers from > 1 to R in a random order. Then you do: > > r = tmp[1] > > So, now r is the first integer from the randomly ordered vector tmp. > Which is the same thing as setting r to be a random integer in 1:R. Then > you do: > > comp_r = setdiff(tmp, r) > > This serts comp_r to be a vector of the integers in tmp, that are not in r. > But since tmp is a permutation of 1:R every element is unique this is the > same as comp_r = tmp[2:R]. Then: > > comp_r = comp_r[1] > > So now comp_r is effectively the second element of the random permutation > of 1:R. What you've effectively done is pick a random integer in the > range 1 to R, but you've gone through a lot of computations to do so. And > you do this every iteration. I think you have made a mistake somewhere, > because you then go on to: > > for k in comp_r > > Since comp_r is a integer this loop will run only once with k taking on > the value of comp_r. The first thing you do in the loop is: > > ind_k = find(comp_r == k) > > which doesn't make any sense since both must be the same value and you > will always get a vector with a single 1 in it. I assume comp_r is > supposed to be a vector of something, but I don't know enough about the > algorithm you are trying to implement to tell what it is your trying to do. > > On Sunday, November 30, 2014 9:29:14 AM UTC+2, Napoleon Vargas wrote: >> >> Dear Julia users, >> >> I'm new to Julia, but from what I've been reading it looks like it has a >> very bright future. I've been working on a function in R that performs >> Metropolis-within Gibbs to estimate linear mixing proportions. The code in >> R works well, but I wanted to speed things up, since I have to try >> different specifications, data, etc., and redoing analysis with several >> data points takes some time. I attached my current R function (along with >> simulated data), and what I think is the "equivalent" Julia code. It is not >> yet complete (it is not wrapped in a function yet), as there are some >> pieces of code that I couldn't quite figure out (I commented out some lines >> in order to make it work temporarily, but the results are not what they're >> supposed to). The R version right now takes around 11 seconds for 25000 >> iterations on a single observation, but for my validation I would probably >> have to try it on several hundred observations, thus my wanting to speed >> things up. >> >> I'd appreciate any help or suggestion. >> >> Napo. >> >
