2
$\begingroup$

I'm trying to understand Bayes Inference better and was wondering if it is possible to use it to fit a function $f(t)$ where I only know the value $f(t)$ for a few values of $t$. For each value of $t$ I have an ensemble of possible values for $f(t)$. As a example I have an experiment which I can repeat as often as I want but my quantaty of interest can only be measured at a few fixed points in time.

As an example lets say $f(t) = sin(\alpha x)$. With $\alpha=1.4$ and the ensemble of each measurment is $Norm(\mu=f(t), \sigma=.01)$. I measure the $f(t)$ at 8 equily distributed points between 0 and $\pi$.

I tried to model this in pymc with:

alpha = pm.Normal('alpha', 3, 1)
tau = pm.Normal('tau', 5, 10)
obs = [pm.CommonDeterministics.Lambda('obs_{}'.format(i), 
                lambda alpha=alpha : np.sin(alpha*xx))
       for i, xx in enumerate(x)]
noises = [pm.Normal('noise_{}'.format(i), o, tau,
                    value=m, observed=True) for i, (m, o) in 
          enumerate(zip(measurements, obs))]

This assumes that the priors for $\alpha$ and $\sigma$ (note: pymc uses $\tau = \frac{1}{\sigma^2}$ instead that is the reason I write $\tau$ in the code.)

When I run this with the MarkovChainMonteCarlo algorithm from pymc I get the following result.

enter image description here

Which is totally off. It does this because it starts to estimate an unreasonably large error $\sigma$ in the observations.

So my questions are.

  1. Is this because I choose bad priors?
  2. Is so can what would be more appropriate priors?
  3. Is something like this in general possible using bayes?

Update

The be clear. I would like to infer the value of $\alpha$ from the observed data.

The whole example can be found here

$\endgroup$
  • $\begingroup$ I think you are estimating independent error terms for each of the measurements. I guess the model you want for your data $y$ is $y \sim N( sin(\alpha x), \tau)$ with your priors $\alpha \sim N(3,1)$, $\tau \sim N(5,10)$ (I don't think these priors are good). So you $\tau$ should be independently drawn from the same distribution for each measurement. $\endgroup$ – thias Jun 5 '15 at 12:21
1
$\begingroup$

Yes, there are plenty of ways of modelling functions. I'd point you in the direction of Gaussian processes for quite a general framework.

Here is an example of modelling a sine wave (green) with one (red) and credible intervals (orange):

enter image description here

Here (IIRC) is the code used to generate it:

from pylab import *
n = 150

ns = 145

X = linspace(-8*pi, 8*pi, n)
S = 2*sin(X)
Y = S + randn(n)
F = linspace(-7*pi,7*pi,ns)

sig1 = 1.0
sig2 = 0.5
b = 0.5

def K(p,q):
  noise = 0
  if p == q:
    noise = sig1
  return sig2*exp(-b*(p-q)**2) + noise

def makeKmn(Xm,Xn):
  m = len(Xm)
  n = len(Xn)
  Kmn = matrix(zeros((m,n)))
  for i in range(n):
    for j in range(n):
      Kmn[i,j] = K(Xm[i],Xn[j])
  return Kmn

Koo = makeKmn(X,X)
Kos = makeKmn(X,F)
Kss = makeKmn(F,F)

YF = Kos.T * Koo.I * matrix(Y).T
YF = array(YF.T)[0]

VF = Kss - Kos.T * Koo.I * Kos
Vs = diag(VF)**(.5)

YFu = 1.96 * Vs + YF
YFl = -1.96 * Vs + YF

fill_between(F,YFu,y2=YFl,facecolor="orange",alpha=0.2)
scatter(X,Y)
plot(F,YF,color="r")
plot(X,S,color='g')
show()

Have a look at gaussianprocess.org for details.

EDIT:

If you suspect function is close to a sine wave one way to proceed is by modelling the mean function of the GP as something like $\mu(x)=b \sin(ax + c)$. If $b$ comes out near zero I expect that would suggest a sine wave is not a good model for your data.

The reason the GP formalism is useful here is that it makes sense to then talk about the value of the function at time points you have never taken a measurement at, i.e. to interpolate. The covariance can be kept simple if desired, e.g. $K(t,t')=\mathbb I_{t=t'}\sigma^2$ (although my hunch would be keeping it so simple isn't going to tell you so much about whether your sine model is good or not.)

$\endgroup$
  • 1
    $\begingroup$ I've read a little bit about gaussian processes now. As far as I understand there are a good tool to make predictions about future observed values based on the available data. While useful this is not what I want to know. I usually know the function my data should follow but not the exact parameters that fit best (the phase $\alpha$ of the sine-curve eg). I don't see how it is possible to infer the value of $\alpha$ with gaussian processes modeling. Are gaussian processes the right tool for that? If so can you point me to some literature. Thanks $\endgroup$ – kain88 Jun 8 '15 at 9:26
  • $\begingroup$ Added some comments. Also, the page linked above has info on software packages for fitting. $\endgroup$ – conjectures Jun 8 '15 at 14:10

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.