Posterior parameter estimation for the Jansen rit

32 views
Skip to first unread message

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.



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()

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


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
 

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 :

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

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


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
 
Reply all
Reply to author
Forward
0 new messages