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

[Tim Leslie]

On 10/26/06, Martin Spacek <scipy at mspacek.mm.st> wrote:
> Hello,
>
> 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 could be wrong here, but doesn't the Ising model only consider
nearest neighbour interactions? Would changing your setup to only
consider these interactions make things work?

Tim

> 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):
> <snip>
>    File
> "C:\bin\Python24\lib\site-packages\scipy\maxentropy\maxentropy.py", line
> 221, in fit
>      callback=callback)
>    File "C:\bin\Python24\Lib\site-packages\scipy\optimize\optimize.py",
> line 811, in fmin_cg
>      callback(xk)
>    File
> "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?
>
> Thanks,
>
> Martin
>
>
> 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:
>
> http://www.nature.com/nature/journal/v440/n7087/abs/nature04701.html
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>