[SciPy-user] Maximum entropy distribution for Ising model - setup?

Martin Spacek scipy at mspacek.mm.st
Thu Oct 26 01:43:24 CDT 2006


This is mostly directed to Ed Schofield, but maybe someone else might be 
able to help as well.

I'm trying to use the scipy.maxentropy model but I've run into a snag, 
and I suspect I simply haven't set it up correctly for my situation.

I'd like to find the maximum entropy distribution of an Ising model of 
magnetic spins (in my case, the spins are actually activity states of 
neurons in a cortical network). Spins can either be down or up (-1 or 1 
respectively). Spins can interact with each other in a pairwise way.

 From what I understand, the Ising model is this:

P(s1, s2, ..., sN) = 1 / Z * exp( sum(hi*si) + 0.5*sum(Jij*si*sj) )
                                    i               i!=j

The si are the spins. I want to model roughly N=10 of them. The hi and 
Jij are the Lagrange multipliers. The hi are the local magnetic fields 
acting on each spin, and the Jij are the degree in which each pair of 
spins interacts. Z is the normalizing partition function.

I'd like to be able to provide the average value for each spin in my 
system (a list of 10 floats in the range (-1, 1)), as well as the 
averages of all possible pairwise products of spins (a list of 10 choose 
2 == 45 floats, also in the range (-1, 1)). Then, after running 
(whichever) maxentropy algorithm and fitting it to the provided data, 
I'd like to be able to pull out the probabilities of each and every 
possible state of the population of spins (say [-1, -1, 1, ...] or [1, 
-1, -1, ...]).

I used scipy/maxentropy/examples/bergerexample.py to get started. I've 
got the following code:

import neuropy as np
from scipy import maxentropy

def f1(x):
     return x==1

class Ising(object):
     """Ising maximum entropy model"""
     def __init__(self, means, pairmeans, algorithm='CG'):
         nbits = len(means)
         samplespace = [-1, 1] # non spiking or spiking
         f1s = [ f1 ]*nbits

         # Return True if they're either both spiking or both non-spiking
         # (correlated), otherwise, return False
         f2s = [ lambda x: f1s[i](x) == f1s[j](x) \
                 for i in range(0, nbits) for j in range(i+1, nbits) ]
         f = np.concatenate((f1s, f2s))
         self.model = maxentropy.model(f, samplespace)
         self.model.verbose = True

         # Now set the desired feature expectations
         means = np.asarray(means)
         pairmeans = np.asarray(pairmeans)
         pairmeans /= 2.0 # add the one half in front of each coefficient
         K = np.concatenate((means, pairmeans))

         # Fit the model
         self.model.fit(K, algorithm=algorithm)

         # Output the distribution
         print "\nFitted model parameters are:\n" \
                + str(self.model.params)
         print "\nFitted distribution is:"
         p = self.model.probdist()
         for j in range(len(self.model.samplespace)):
             x = self.model.samplespace[j]
             print '\tx:%s, p(x):%s' % (x, p[j])

i = Ising(means, pairmeans)

means is a list of 10 floats that mostly hover around -0.8 to -0.5. 
pairmeans is a list of 45 floats that are mostly around 0.4 to 0.9. This 
runs, but gives me a DivergenceError:

Grad eval #0
   norm of gradient = 6.02170024163
Function eval # 0
   dual is  0.69314718056
Function eval # 1
   dual is  -29.5692018412
Grad eval #1
   norm of gradient = 5.03761439516
Function eval # 2
   dual is  -147.846016909
Grad eval #2
   norm of gradient = 5.03761183138
Function eval # 3
   dual is  -620.953271018
Grad eval #3
   norm of gradient = 5.03761183138
Function eval # 4
   dual is  -1478.46016909
Grad eval #4
   norm of gradient = 5.03761183138
Iteration # 0
Function eval # 5
   dual is  -1478.46016909
Traceback (most recent call last):
"C:\bin\Python24\lib\site-packages\scipy\maxentropy\maxentropy.py", line 
221, in fit
   File "C:\bin\Python24\Lib\site-packages\scipy\optimize\optimize.py", 
line 811, in fmin_cg
"C:\bin\Python24\lib\site-packages\scipy\maxentropy\maxentropy.py", line 
397, in log
     raise DivergenceError, "dual is below the threshold 'mindual'" \
DivergenceError: "dual is below the threshold 'mindual' and may be 
diverging to -inf.  Fix the constraints or lower the threshold!"

I tried changing the model's mindual attrib from -100 to something 
ridiculous like -100000, but that didn't change much, just took more 
iterations before bombing out. I've also tried some of the algorithms 
other than 'CG', without any luck.

I don't know much about maximum entropy modelling, especially the actual 
algorithms. I only have a general idea of what it's supposed to do, so 
the above code is really just a shot in the dark. I suspect the problem 
is in properly defining both the means and the pairwise means for the 
model. Or it may have something to do with continuous vs. discrete 
values. Any suggestions?



P.S. For anyone interested, this is related to the recent Nature paper 
"Weak pairwise correlation imply strongly correlated network states in a 
neural population" by Schneidman et al:


More information about the SciPy-user mailing list