Hi, I am trying to make a bayesian fit of a microlensing light curve, I prepare these model and i adapted some parts using jax:
def model_d(observado,tiempo):
w_min=np.where(np.min(observado)==observado)[0][0]
t_0=numpyro.sample("t_0", dist.Uniform(tiempo[w_min]-500,tiempo[w_min]+500))
u_0=numpyro.sample("u_0", dist.Uniform(1e-30,1.0))
t_E=numpyro.sample("t_E", dist.Uniform(100.,300.0))
q=numpyro.sample("q", dist.Uniform(1e-30,1.0))
s=numpyro.sample("s", dist.Uniform(1e-30,2))
alpha=numpyro.sample("alpha", dist.Uniform(1e-30,2*jnp.pi))
tau = (tiempo-t_0)/t_E
beta = u_0
x_old,y_old,theta=tau,beta,alpha
#binrot(alpha, tau, beta)
cos_theta = jnp.cos(theta)
sin_theta = jnp.sin(theta)
x = x_old * cos_theta - y_old * sin_theta
y = x_old * sin_theta + y_old * cos_theta
x = x - s*q/(1.+q)
zeta = x + 1j*y
cofs0=(1+q)**2*(s+zeta.conjugate())*zeta.conjugate()
cofs1=(1+q)*(s*(q-abs(zeta)**2*(1+q))+(1+q)*((1+2*s**2)-abs(zeta)**2+2*s*zeta.conjugate())*zeta.conjugate())
cofs2=(1+q)*(s**2*q-s*(1+q)*zeta+(2*s+s**3*(1+q)+s**2*(1+q)*zeta.conjugate())*zeta.conjugate()-2*abs(zeta)**2*(1+q)*(1+s**2+s*zeta.conjugate()))
cofs3=-(1+q)*(s*q+s**2*(q-1)*zeta.conjugate()+(1+q+s**2*(2+q))*zeta+abs(zeta)**2*(2*s*(2+q)+s**2*(1+q)*(s+zeta.conjugate())))
cofs4=-s*(1+q)*((2+s**2)*zeta+2*s*abs(zeta)**2)-s**2*q
cofs5=-s**2*zeta
coefs = jnp.vstack((cofs0, cofs1,cofs2,cofs3,cofs4,cofs5))
z0=jnp.apply_along_axis(jnp.roots,0,coefs,strip_zeros=False)
W1 = (1./(1.+q)*(1./z0+q/(z0+s)))
z1=z0.copy()
z1 =jnp.where(jnp.abs(z0-W1.conjugate()-zeta)>0.01, jnp.nan,z1)
W2 = -1./(1.+q)*(1./z1**2+q/(z1+s)**2)
W2=jnp.nan_to_num(W2,1e-30)
mu0 = 1./(1.-jnp.abs(W2)**2)
mu1=mu0.copy()
A = jnp.sum(jnp.abs(mu1),axis=0)
mag0=numpyro.sample("mag0", dist.Uniform(17,23))
fs=numpyro.sample("fs", dist.Uniform(1e-30,1.0))
mu=mag0-2.5*jnp.log10(fs*(A-1)+1)
numpyro.sample("mag", dist.Normal(mu, 0.001), obs=observado)
How you can observe in some parts i had to fine the roots of a polynom, whit complex numbers, then to run the fit I use:
#Running a HMC-NUTS
with numpyro.validation_enabled():
rng_key = random.PRNGKey(20)
rng_key, rng_key_ = random.split(rng_key)
num_warmup, num_samples = 4000, 10000
kernel = NUTS(model_d, init_strategy=init_to_median(), target_accept_prob=0.9)
mcmc = MCMC(kernel, num_warmup=num_warmup, num_samples=num_samples)
mcmc.run(rng_key_,observado=observado,tiempo=tiempo)
mcmc.print_summary().
After run, this error appears “Cannot find valid initial parameters. Please check your model again.”
I tried to look for answers in the forum, but i don’t find nothing similar and also i dint know if the code is the problem or if the problem is a stats miss conception of my self, so I really appreciate a little light on these.
Thanks Felipe
PD: observados= an array of values from 18-23 and tiempo: another array this time whit values in the order of 2458762, and the version of numpyro is 0.10.