![](https://secure.gravatar.com/avatar/ac4e2152839c56c5326cd95ae54296e9.jpg?s=120&d=mm&r=g)
On Thu, Jan 7, 2010 at 9:19 AM, Jordi Molins Coronado <jordi_molins@hotmail.com> wrote:
Sorry, I see my previous message has been a disaster in formatting. I try now in a different way. Sorry for the inconveniences.
Hello, I am new to this forum. I am looking for a numerical solution to the inverse problem of an Ising model (or a model not-unlike the Ising model, see below). I have seen an old discussion, but very interesting, about this subject on this forum (http://mail.scipy.org/pipermail/scipy-user/2006-October/009703.html). I would like to pose my problem (which is quite similar to the problem discussed in the thread above) and kindly ask you your opinion on that:
My space is a set of discrete nodes,s_i, where i=1,...,N, which can take two values, {0,1}. Empirically I have the following information: <s_i>_emp and <s_i*s_j>_emp, where i,j=1,...,N with i!=j.
It is well known in the literature that the Ising model
P(s_1, s_2, ..., s_N) = 1 / Z * exp( sum{for all i}(h_i*s_i) + 0.5*sum{for all i!=j}(J_ij*s_i*s_j) ) maximizes entropy with the constraints given above (in fact, this is not the Ising model, because the Ising model assumes only nearest-neigbour interactions, and I have interactions with all other nodes, but I believe it is still true that the above P(s1,...sN) still maximizes entropy given the constraints above). What I would like is to solve the inverse problem of finding the h_i and J_ij which maximize entropy given my constraints. However, I would like to restrict the number of h_i and J_ij possible, since having complete freedom could become an unwieldly problem. For example, I could restrict h_i = H and J_ij = J for all i,j=1,...N, i!=j, or I could have a partition of my nodes, say nodes from 1 to M having h_i = H1 and J_ij=J1 i,j=1,...,M i!=j, and h_i=H2 and J_ij=J2 i,j=M+1,...,N i!=j, and the J_ij=J3 for i=1,...,M and j=M+1,...N. If I understand correctly the discussion in the thread shown above, a numerical solution for the inverse problem would be: hi_{new}=hi_{old} + K * (<si> - <si>_{emp}) Jij_{new}=Jij_{old}+ K' * (<si*sj> - <si*sj>_{emp})
where K and K' are pos. "step size" constants. (On the RHS, <si> and <si*sj> are w.r.t. hi_{old} and Jij_{old}.) Have I understood all this correctly? In particular, for the case h_i = H and J_ij = J for all i,j=1,...N, i!=j could I simplify the previous algorithm by restricting the calculations only to say i=1 (i=2,...,N should be the same?), and for the case h_i = H1 and J_ij=J1 i,j=1,...,M i!=j, and h_i=H2 and J_ij=J2 i,j=M+1,...,N i!=j simplify it by restricting the calculations only to say i=1 and i=M+1? Thank you for your help and sorry if I am new here and I have committed some "ettiquette" mistake.
Hi, I'm not so familiar with the statisitical mechanics notation, but you might be interested in the maxent module of the pyentropy package I have produced as part of my PhD: http://code.google.com/p/pyentropy/ The main purpose of pyentropy is calculation of bias corrected entropy and information values from limited data sets, but it includes the maxent module which computes maximum entropy distributions over finite alphabet spaces with marginal constraints of up to any order. (I am working in computational neuroscience so much of the notation will probably be a bit different). So I think a second order solution from this framework over a binary space is the same as the Ising model. You can get the hi and J's directly from the results (they are called theta in the code) - although I think they have a slightly different normalisation because of Ising being -1,1 and this being 0,1... http://pyentropy.googlecode.com/svn/docs/api.html#module-pyentropy.maxent With this on a normal computer I can solve for about 18 binary variables in a reasonable amount of time (ie less than an hour for several runs), but it becomes highly exponential with more vectors. (I have a much more efficient but more hackish version of the same algorithm that I haven't added to pyentropy yet, but will in the next few weeks). In the case you describe where the thetas are constrained to be equal at each order the system can be solved much more efficiently. I have code to do this which is not released (and a bit messy) but if you are interested I could send it to you. In neuroscience this situation is called the 'pooled model'. You can see a description for how to solve in this reduced case here: http://rsta.royalsocietypublishing.org/content/367/1901/3297.short The method is using information geomery from Amari - by transforming between the P space (probability vector), eta space (marginals) and theta space (the h,J's etc.) it is possible to find the maximum entropy solution as a projection. Anyway I'm not sure if that helps... probably the documentation on how to get the thetas and the ordering of the vector might not be so clear, so give me a shout if you start using it and have any questions. If you don't have the full probability distribution you can pass the eta vector to the solve function (which would be a vector of <si>, <si,sj> ) and set optional argument eta_given=True (hopefully clear from the option handling in the code). Cheers Robin