Posterior parameter estimation for the Jansen rit

Deepa Tilwani

unread,

Nov 21, 2022, 2:54:10 PM11/21/22

to TVB Users

how to define priors before the jansen rit, and simulate from the jansen rit and then estimate the posterior parameters using MCMC ?
Basically inversion on simulated sEEG/EEG parameters.

WOODMAN Michael

unread,

Nov 22, 2022, 3:09:13 AM11/22/22

to TVB Users

hi,

IIRC your previous email was about the forward model (which for s/EEG is a gain matrix that should be calculated before hand with OpenMEEG or Brainstorm or similar software).

The particular problem of inverting EEG with the Jansen Rit has been well explored within the DCM framework and it's probably a good idea to understand the content of this paper before proceeding

https://www.sciencedirect.com/science/article/pii/S1053811921009356

cheers,

Marmaduke

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Deepa Tilwani <tilwan...@gmail.com>
Sent: Monday, November 21, 2022 8:54:10 PM
To: TVB Users
Subject: [RESEAUX SOCIAUX] [TVB] Posterior parameter estimation for the Jansen rit

how to define priors before the jansen rit, and simulate from the jansen rit and then estimate the posterior parameters using MCMC ?
Basically inversion on simulated sEEG/EEG parameters.

--
You received this message because you are subscribed to the Google Groups "TVB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tvb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/c64ac929-f087-4ea0-b795-a7abbb40749bn%40googlegroups.com.

Deepa Tilwani

unread,

Nov 26, 2022, 1:09:25 PM11/26/22

to tvb-...@googlegroups.com

Hi Thanks for your email.

I was trying to customize the parameters of the Jansen rit model by doing

#observations model to call in the mcmc

def jrm_model(A, B, a, b , J, a_1, a_2, a_3, a_4, vmax, v0):

#Calling inbuilt model
jrm = models.JansenRit()

jrm.A = jnp.asarray(A)

jrm.B = jnp.array([B])

jrm.a =jnp.array([a])

jrm.a_1 =jnp.array([a_1])

jrm.a_2 =jnp.array([a_2])

jrm.a_3 =jnp.array([a_3])

jrm.a_4 =jnp.array([a_4])

jrm.b =jnp.array([b])

jrm.J =jnp.array([J])

jrm.v0 =jnp.array([v0])

jrm.nu_max = jnp.array([vmax])

phi_n_scaling = (jrm.a * jrm.A * (jrm.p_max-jrm.p_min) * 0.5 )**2 / 2.
sigma = np.zeros(6)
sigma[3] = phi_n_scaling

# the other aspects of the simulator are standard

# Sim Initialise Simulator.
sim = simulator.Simulator(
model=jrm, # Initialise the Jansen rit model
connectivity=connectivity.Connectivity.from_file(),# Initialise a Connectivity.
coupling=coupling.SigmoidalJansenRit(),
integrator=integrators.HeunStochastic(dt=2 ** -4, noise=noise.Additive(nsig=sigma)),
monitors=(monitors.TemporalAverage(period=2 ** -2),),
simulation_length=1e3,
).configure()

# run it -- genrate data
(times ,seeg), = sim.run()
print(seeg.shape)

return seeg

But it's is giving me error :

Attribute can't be set to an instance of <class 'jaxlib.xla_extension.DeviceArray'>
  attribute tvb.simulator.models.jansen_rit.JansenRit.A = NArray(label=':math:`A`', dtype=float64, default=array([3.25]), dim_names=(), ndim=None, required=True)

Can someone help me out if they know so. I tried changing the dtype but it is not working.

If I am calling inbuilt function without customizing any parameters it's working fine.

Deepa Tilwani,

PhD Student,

UofSC Email : dtil...@mailbox.sc.edu

Laboratory for Integrative Neuroscience Analysis,

Artificial Intelligence Institute (AIISC),

Carolina Autism and Neurodevelopment (CAN) Research Center,

Institute for Mind and Brain (IMB)

LinkedIn | https://deepa-tilwani.github.io/

You received this message because you are subscribed to a topic in the Google Groups "TVB Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/tvb-users/pKjMD9H8w1s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to tvb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/5e85086c121a4088afdd34fa7698b7ef%40univ-amu.fr.

Deepa Tilwani

unread,

Nov 27, 2022, 10:53:32 PM11/27/22

to tvb-...@googlegroups.com

I tried passing in jrm.A = np.array([A]) It is giving me this error.


The numpy.ndarray conversion method __array__() was called on the JAX Tracer object Traced<ConcreteArray(-0.8757257461547852, dtype=float32)>with<JVPTrace(level=2/0)> with
  primal = DeviceArray(-0.87572575, dtype=float32)
  tangent = Traced<ShapedArray(float32[])>with<JaxprTrace(level=1/0)> with
    pval = (ShapedArray(float32[]), None)
    recipe = LambdaBinding()

Deepa Tilwani,

PhD Student,

UofSC Email : dtil...@mailbox.sc.edu

Laboratory for Integrative Neuroscience Analysis,

Artificial Intelligence Institute (AIISC),

Carolina Autism and Neurodevelopment (CAN) Research Center,

Institute for Mind and Brain (IMB)

LinkedIn | https://deepa-tilwani.github.io/

WOODMAN Michael

unread,

Nov 28, 2022, 3:04:31 AM11/28/22

to tvb-...@googlegroups.com

Hi,

First, thanks for making an effort with TVB, the API is understandably not what you'd expect if you're accustomed to Jax or Tensorflow, and it does some type checking, which you encountered before. Since our full-feature simulator uses some NumPy APIs which won't be compatible with tracer-based systems like Jax, you will need to use the code generation backend to obtain some plain NumPy code, which can then be traced with Jax. There's an effort to do this also for Theano (Aesara from PyMC devs), going on in a PR here

https://github.com/the-virtual-brain/tvb-root/pull/630

Something similar could be done for Jax, with some tweaks to some of the API calls. But, it also seems that aesara supports Jax as a backend, so you may just try that branch directly.

Lastly, the Jansen-Rit model is not (yet) available for code generation, but it is just missing some metadata on the model class. You can see here

https://github.com/the-virtual-brain/tvb-root/pull/630/commits/aebfc85d61fd39d099ac05118c66d48ed084b632#diff-be317ee6b12688f97cf40d9d34d6f5769441c3720c3e8a2be73e51cd47d56a6eR331

what metadata is required, and it would be more than welcome to have that contribution for the Jansen Rit model.

cheers,

Marmaduke

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Deepa Tilwani <tilwan...@gmail.com>

Sent: Monday, November 28, 2022 4:52:54 AM
To: tvb-...@googlegroups.com
Subject: [RESEAUX SOCIAUX] Re: [RESEAUX SOCIAUX] [TVB] Posterior parameter estimation for the Jansen rit

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/CAFHBvivAwiu3Y6UncPMjG4y4UTNwnKc7Ysa4M%3D5isV%2BonD-dLQ%40mail.gmail.com.

Deepa Tilwani

unread,

Nov 28, 2022, 6:20:31 PM11/28/22

to tvb-...@googlegroups.com

Hi,

Thanks for the earlier email.

I changed a bit for passing values, here jnp is imported as: import jax.numpy as jnp

def jrm_model(A, B, a, b , J, a_1, a_2, a_3, a_4, vmax, v0):

#Calling inbuilt model

jrm = models.JansenRit()

print(np.array([jnp.ravel(A)]))
print(type(np.array([jnp.ravel(A)])[0]))

jrm.A = np.array([jnp.ravel(A)])

"""

jrm.B = jnp.array([B])

jrm.a =jnp.array([a])

jrm.a_1 =jnp.array([a_1])

jrm.a_2 =jnp.array([a_2])

jrm.a_3 =jnp.array([a_3])

jrm.a_4 =jnp.array([a_4])

jrm.b =jnp.array([b])

jrm.J =jnp.array([J])

jrm.v0 =jnp.array([v0])

jrm.nu_max = jnp.array([vmax])

"""

phi_n_scaling = (jrm.a * jrm.A * (jrm.p_max-jrm.p_min) * 0.5 )**2 / 2.
sigma = np.zeros(6)
sigma[3] = phi_n_scaling

#print(jrm)

# the other aspects of the simulator are standard

# Sim Initialise Simulator.
sim = simulator.Simulator(
model=jrm, # Initialise the Jansen rit model
connectivity=connectivity.Connectivity.from_file(),# Initialise a Connectivity.
coupling=coupling.SigmoidalJansenRit(),
integrator=integrators.HeunStochastic(dt=2 ** -4, noise=noise.Additive(nsig=sigma)),
monitors=(monitors.TemporalAverage(period=2 ** -2),),
simulation_length=1e3,
).configure()

# run it -- genrate data
(times ,seeg), = sim.run()
print(seeg.shape)

return seeg

for calling : jrm_model(0, 1,1,1,1,1,1,1,1,1,1)

I encountered with another error :

If someone knows why it is giving shapes error, please let me know;

Deepa Tilwani,

PhD Student,

UofSC Email : dtil...@mailbox.sc.edu

Laboratory for Integrative Neuroscience Analysis,

Artificial Intelligence Institute (AIISC),

Carolina Autism and Neurodevelopment (CAN) Research Center,

Institute for Mind and Brain (IMB)

LinkedIn | https://deepa-tilwani.github.io/

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/242b19fbe7d34031beacd28d9ee68298%40univ-amu.fr.

WOODMAN Michael

unread,

Nov 29, 2022, 3:32:10 AM11/29/22

to tvb-...@googlegroups.com

hi

TVB uses a last extra dimension for the modal models which JR is not. You can just add the extra dimension to your array so that the shape matches.

cheers,

Marmaduke

From: tvb-...@googlegroups.com <tvb-...@googlegroups.com> on behalf of Deepa Tilwani <tilwan...@gmail.com>

Sent: Tuesday, November 29, 2022 12:19:53 AM
To: tvb-...@googlegroups.com
Subject: [RESEAUX SOCIAUX] Re: [RESEAUX SOCIAUX] Re: [RESEAUX SOCIAUX] [TVB] Posterior parameter estimation for the Jansen rit

To view this discussion on the web visit https://groups.google.com/d/msgid/tvb-users/CAFHBviso9pbROKFoy46XMbML3mXPB4tSZL3PgJX1ebZq3cOU%2Bg%40mail.gmail.com.

Reply all

Reply to author

Forward