[SciPy-user] Maximum entropy distribution for Ising model - setup?
Martin Spacek
scipy at mspacek.mm.st
Thu Oct 26 01:43:24 CDT 2006
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'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
More information about the SciPy-user
mailing list