[SciPy-user] Maximum entropy distribution for Ising model - setup?
Tim Leslie
tim.leslie at gmail.com
Fri Oct 27 00:45:38 CDT 2006
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
>
More information about the SciPy-user
mailing list