1
$\begingroup$

I've asked this question on stackoverflow too, but no answer yet. This seems a more appropariate place to ask this question:

I'm messing around with pymc to understand it a bit better. Now I am trying to fit a logistic function to my generated data. Though the parameter it returns is not even close to the real parameter. Any idea why ?

import pymc as pm
from matplotlib import pyplot as plt
import numpy as np
from pymc.Matplot import plot as mcplot


#data generation
def logistic(x, beta):
    return 1.0 / (1.0 + np.exp(beta * x))

#hidden
beta = -5

x = np.linspace(-4, 4, 600)
data = logistic(x, -5)

plt.plot(x,data , label=r"$\beta = -5$")
plt.legend();

plt.show()




#Priors
beta = pm.Uniform("beta", -10, 10, value=1.)
sig = pm.Uniform("sig", 0.0, 100.0, value=1.)


@pm.deterministic
def logistic(x=data, beta=beta):
    return 1.0 / (1. + np.exp(beta * x))


#likelihood
y = pm.Normal("obs", mu=logistic, tau=1.0/sig**2, value=data, observed=True)


model = pm.Model([y, data,beta,sig])
mcmc = pm.MCMC(model) 
mcmc.sample(20000, 5000, 1)


beta_samples = mcmc.trace('beta')[:]
sig_samples = mcmc.trace('sig')[:]

print beta_samples.mean()
print sig_samples.mean()

mcplot(mcmc.trace('beta'))
mcplot(mcmc.trace('sig'))
plt.show()
$\endgroup$
  • $\begingroup$ As the help explains, please don't crosspost. Choose one site that's best suited to your question and post it there. If you later change your mind you can flag it and ask for it to be moved. $\endgroup$ – Glen_b -Reinstate Monica Apr 1 '14 at 11:44
1
$\begingroup$

Your error is in:

@pm.deterministic
def logistic(x=data, beta=beta):
    return 1.0 / (1. + np.exp(beta * x))

You are passing in x=data, but data is the sigmoid'ed from above. Try this:

@pm.deterministic
def logistic(x=x, beta=beta):
    return 1.0 / (1. + np.exp(beta * x))

That should work.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.